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L  INTRODUCTION 


A.  PROBLEM  STATEMENT 

The  dynamic  tesponse  of  a  submeraible  vehicle  operating  at  the  extremes 
of  its  operational  envelope  is  becoming  increasingly  important  in  order  to  en¬ 
hance  vehicle  operations.  Traditionally,  (fynamic  stability  of  motion  is  stud¬ 
ied  using  eigenvalue  analysis  where  the  equations  of  motion  are  linearised 
around  nominal  straight  line  level  ffight  paths  (Arentsen  &  Mandel,  1960), 
(Clayton  &  Bishop,  1982),  (Feldman,  1987).  Under  certain  simplified  as¬ 
sumptions,  a  simple  criterion  >  0  can  be  obtained  where  the  stability 
index  is  function  of  the  hydrodynamic  coefficients  in  heave  and  latch. 
Values  for  the  stability  index  can  be  computed  by, 

^  ,  Jll.(Z, -l-m) 

~z,Mr  ■ 

This  index  is  analogous  to  the  familiar  stability  coefficient  for  horisontal 
plane  maneuvering  (Lewis,  1989)  and  can  be  thought  of  as  a  high  speed  iq>- 
proximation  where  the  effect  of  the  metacentric  restoiizig  moment  is  minimal. 
If  the  value  of  6r«  is  greater  than  zero,  the  vehicle  is  dynamically  stable.  As 
we  point  out  in  the  next  chapter  though,  this  is  only  a  sufficient,  and  rather 
conservative  condition  for  stability.  It  is  not  a  necessary  condition  in  the 
sense  that  <  0  indicates  instability  at  infinite  forward  speed.  It  is  quite 
possible  that  at  normal  <^erating  speeds  the  vehicle  might  be  directionally 
stable.  Furthermore,  <  0  indicates  a  divergent  loss  of  stability  which  is 
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quite  uncommon  in  the  vertical  plane.  Moat  modem  aubmarinea  exhibit  a 
flutter-like  inatability  at  high  apeed,  which  can  not  be  analyzed  uaing  the 
above  aimplifled  index.  Divergent  motiona  may  develop  in  combined  aix 
degreea  of  freedom  (Papouliaa  et  al  1993)  and  their  occurrence  can  not  be 
analyzed  by  a  aingle  at  ability  index. 


B.  OBJECTIVES  AND  OUTLINE 

In  thia  work  we  examine  the  problem  of  at  ability  of  motion  with  controla 
fixed  in  the  vertical  jdane,  with  particular  emphaaia  on  the  mechaniam  of 
loaa  of  atability  atraight  line  motion.  The  aurge  equation  ia  decoupled 
from  heave/pitch  through  a  perturbation  aeriea  approach  (Bender  &  Orazag, 
1978).  It  ia  ahown  that  loaa  of  atability  occura  in  the  form  of  generic  bifurca- 
tiona  to  periodic  solutiona  (Guckenheimer  &  Holmea,  1983).  Ihylor  expan* 
siona  and  center  manifold  approximationa  are  employed  in  order  to  iaolate 
the  main  nonlinear  terma  that  influence  ayatem  reaponae  after  the  initial  loaa 
atability  (Haaaard  &  Wan,  1978).  Integral  avera^ng  ia  performed  in  order 
to  combine  the  nonlinear  terma  into  a  deaign  atability  coefficient  (Chow  & 
Mallet-Paret,  1977).  Special  attention  ia  paid  to  the  atudy  of  the  quadratic 
dr^  terma  aa  they  conatitute  aome  of  the  main  nonlinear  terma  of  the  equa- 
tiona  of  motion.  Thia  ia  a  unique  feature  of  the  open  loop  problem.  In  aimilar 
atudiea  of  loaa  of  atability  under  dewed  loop  depth  control  (Bateman,  1993) 
it  waa  found  that  the  main  damping  mechaniam  ia  provided  through  the  uae 
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of  control  susrafces.  The  difficulty  associated  with  the  nonsmoothness  the 
absolute  value  nonlinearities  of  the  quadratic  drag  forces  is  dealt  with  by 
employing  the  concept  of  generalized  gradient  (Clarke,  1983).  This  has  the 
advantage  of  keeping  the  linear  terms  constant,  unlike  the  linear/cubic  ap¬ 
proximation  typically  used  in  ship  toll  nx>tion  studies  (Dalzell,  1978),  where 
the  linear  damping  coefficient  is  a  function  of  the  assumed  amplitude  of  mo¬ 
tion. 

Vehicle  modeling  in  this  work  follows  standard  notation  (Gertler  &  Hagen, 
1976),  (Smith  et  al  1978),  and  numerical  results  are  presented  for  the  DARPA 
SUBOFF  model  (Roddy,  1990)  for  which  a  set  of  hydrodynamic  coefficients 
and  geometric  properties  is  available.  Furthermore,  the  baseline  vehicle  is 
marginally  stable  with  controls  fixed  under  normal  operating  conditions  and 
can  serve  as  a  prime  example  for  the  techniques  described  in  this  work.  The 
model  has  been  experimentally  validated  for  angles  of  attack  on  the  hull 
between  ±15  deg.,  while  the  constant  coefficient  iq>proximation  introduces 
very  little  error  in  time  domain  amulations  (Tinker,  1978).  Unless  otherwise 
mentioned,  all  results  in  this  work  are  presented  in  standard  dimensionless 
form  with  respect  to  the  vehicle  length  I  =  4.26  m,  and  nominal  forward 
speed  U  =  2.44  m/sec.  All  angular  deflections  are  shown  in  degrees. 
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n.  PROBLEM  FORMULATION 


A.  EQUATIONS  OP  MOTION 

Assuinixig  that  vehicle  motion  is  restricted  in  the  vertical  plane,  the  math¬ 
ematical  model  consists  of  the  coupled  nonlinear  heave  and  pitch  equations 
of  motion.  In  a  moving  coordinate  frame  fixed  at  the  vehicle’s  geometrical 
center,  Newton’s  equations  of  motion  for  a  port /starboard  qrmmetric  and 
neutrally  buoyant  vehicle  are  expressed  in  dimensionless  form  as  follows, 

m(w  -uq  —  za^  —  zaq)  =  Z^q  -f-  Z^w  +  Z^q  +  Z^„w 

/noM 

-Cd  /  b{x)(w  —  sq)|tt>  —  *q|  d*  -f  ZsS  ,  (2) 

Jtall 

lyi  +  mza{u  +  vjq)  -  masaiii}  —  uq)  =  M^q  +  M^w  -t-  M,q  -f  M^w 

+Cd  I  b{x){w  —  xq)\vo  —  xq\x  dx 

ytau 

-xqbW cos —  zqbW Aa0  +  MfS  ,  (3) 

where  xqb  =  —  —  xa,  and  the  rest  the  symbols  are  based 

on  standard  notation  and  are  explained  in  Table  1.  Without  loss  of  generality 
we  can  assume  that  za  =  sea  =  0,  so  that  xqb  =  and  sera  =  The  cross 
flow  integral  terms  in  these  equations  become  very  important  for  high  angles 
of  attack  maneuvering,  where  they  provide  the  primary  motion  damping. 
The  drag  coefficient,  Ca*  is  assumed  to  be  constant  throughout  the  vehicle 
length  for  simplicity.  This  does  not  affect  the  qualitative  properties  of  the 
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results  that  fdlow.  The  vehicle  pitch  rate  is, 

0  =  q.  (4) 

D3maimc  coupling  between  surge  and  heave/pitch  is  present  due  to  coordi¬ 
nate  coupling  as  a  result  of  the  nonzero  metacentric  height.  Therefore,  pitch 
and  heave  motions  must  be  studied  together  with  surge, 

mti  +  imoq  —  mzaq^  -I-  mzaq  =  -h  XiiU  +  X^wq  -h  X^wto^ 
-{■X^uU^  +  X„nn^  +  Xu6\  (5) 

where  we  assume  that  both  resistance  and  propulsive  forces  are  proportional 
to  the  square  of  the  speed  or  the  propeller  revolutions,  respectively. 

In  analyzing  controls  fixed  stability  of  nation,  the  case  ^  =  0  is  examined 
first.  The  steady  state  solutions  of  the  equations  are  determined  by  to  =  g  = 
u  =  0  =  9)  =  0,  where  subscript  0  indicates  variable  value  at  steady  state. 
Substituting  these  conditions  in  (2)  we  get, 

Z«too  —  CD>l«u>o|tool  =  0  ,  (6) 

where, 

A,=/76(*)<te,  (7) 

is  the  ‘Vaterplane”  area.  Since  Z„  <  0,  equation  (7)  admits  only  one  solu¬ 
tion,  namely  Wo  =  0.  Equation  (3)  then  yields, 

tan^o  =  -— ,  (8) 

ZqB 

while  (5)  is  used  to  determine  the  nominal  forward  speed,  uq. 
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Zi,  Z] 
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TABLE  1:  NOMENCLATURE 


dununy  independent  variable 


steady  state  value  of  a 


expansion  coefficients  in  terms  of  si ,  S3 


local  beam  of  the  hull 


quadratic  drag  coefficient 


regularization  parauneter 


stem  {dane  deflection 


perturbation  parameter,  e  = 


criticality  difference,  6  =  u  —  Ue 


vehicle  mass  moment  of  inertia 


nonlinear  stability  coefficient 


system  eigenvalue 


vehicle  mass 


pitch  moment 


derivative  of  M  with  respect  to  a 


pitch  rate  _ 


polar  coordinates  of  Si,  zj 


transformation  matrix  of  x  to  a 


pitch  angle 


forward  speed 


critical  value  of  u 


heave  velocity 


state  variables  vector,  x  = 


surge  force 


derivative  of  X  with  respect  to  a 


body  fixed  coordinates  of  vehicle  center  of  buoyancy 


body  fixed  coordinates  of  vehicle  center  of  gravity 


center  of  gravity/center  of  buoyancy  separation,  xq  —  xg 


vehicle  metacentric  height,  zq  —  zg 


state  variables  vector  in  its  rxirmal  form 


critical  coordinates  of  % 


stable  coordinate  of  z 


heave  force 


derivative  of  Z  with  respect  to  a 


imaginary  part  of  critical  pair  of  eigenvalues 
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B.  REDUCTION  OF  ORDER 

The  linearized  surge,  heave,  and  pitch  equations  o[  motion  in  the  vicinity 
of  the  nominal  point  are, 

(m  —  +  mzGq  =  ,  (9) 

(m  -  Z^)w  -  (mxa  +  Z4)  =  (Z,  +  m)q  +  Z^w  ,  (10) 

(/y  —  M^)q  +  mz<3U  -  (m*<5  +  M^)w  =  M„w  +  (M,  -  mxG)q 

+Mte ,  (11) 

where, 

M$  =  xgb  sin  do  —  W  cos  6q  ,  (12) 

is  the  hydrostatic  restoring  naoment  coefficient.  The  characteristic  equation 
of  (4),  (9),  (10),  and  (11)  is  obtained  as, 

(— i4i A  +  Ri)(i42A*  +  RjA*  +  CjA  +  Dj)  +  i48A^  +  R3A®  =  0  ,  (13) 

where, 

—  tn  ACii  , 

Bx  —  , 

At  =  (m  -  Z^){Iy  -  M^)  -  (Z4  +  mxG)(M^  +  rrixc) , 

Bt  =  -Z„{Iy  -  Mi)  - {m-  Z,b)(Af,  -  itubg)  -  (^g  +  Tn)(Afw  +  mxa) 
-M„{Zi  +mxG) , 

Ct  =  -M${m  -  Z*)  +  Zy,{M^  —  ttixg)  -  My,(Zq  +  m)  , 

Bt  =  MfZyg  , 
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Aa  =  (m*G)’(m  -  Z«)  , 


Ba  =  -{mzaYZy,  . 


It  can  be  seen  that  the  parameter  (mza)  u  responsible  for  surge  and  heave, 
pitch  coupling.  For  zq  =  0,  equation  (13)  decouples  into  the  surge  eigenvalue 
A  =  Bi/Ai  and  the  classical  cubic  characteristic  equation  for  the  vertical 
plane.  It  should  be  mentioned  that  the  effect  of  the  forward  speed  tt  is 
embedded  into  the  definition  for  the  dimensionless  vehicle  weight  W  through, 


W — . 


W 


(M) 


If  we  introduce  a  smallness  parameter. 


e  =  (n«o)’  , 


(15) 


we  can  rewrite  (13)  in  the  form, 

(A  +  ea)\*  +  iB  +  e^)A’  +  CA*  +  27A  +  £?  =  0  ,  (16) 


where  in  terms  of  previously  defined  coefficients, 


A  =  -Ai<h  1 
B  =  —  Ba  -i-  BiAa  , 
C  =  —A\Ca  +  BiBa  , 
D  =  —AxDa  BiCa  « 
E  =  B1D2  y 

a  =  m  —  y 

=  -Z^. 
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Following  (Bender  &  Qrszag,  1978)  we  expand  the  w^utions  of  (16)  in  a 
regular  perturbation  series, 


A  =  Ao  +  XiC  +  , 


(17) 


where  Ao  is  an  eigenvalue  for  e  =  0  (uncoupled  surge  or  heave/pitch),  Ai  is 
the  first  order  correction  due  to  dynamic  coupling,  and  0(e’)  contains  second 
and  higher  order  terms  in  e.  Substituting  (17)  into  (16)  we  can  get, 


X _ Ag(Q!Ao  +fi)  j 

^  ^  4AA2  +  3RA2  4- 2<7Ju -1- D  ^  ^  ' 


It  can  be  seen  that  the  correction  term  is  very  small  compared  to  the  un¬ 
coupled  root  as  evidenced  by  the  Aq  term.  This  is  particularly  true  when 
Ao  nears  zero;  i.e.,  dose  to  a  bifurcation  point.  Therefore,  loss  of  stability 
can  be  studied  analyzing  the  heave/pitch  equations  decoupled  from  surge. 
The  characteristic  equation  then  becomes, 


i4jA*  +  RjA*  -f-  C^jA  -f-  =  0  . 


(19) 


Plots  of  the  system  eigenvalues  at  nominal  speed  versus  zq  are  shown  in 
Figures  1  through  4.  The  surge  eigenvalues  is  real  and  negative  throughout 
the  range  of  zot  while  the  heave/pitch  eigenvalues  are  real  for  small  values 
of  zq.  The  two  larger  real  heave/pitch  eigenvalues  coalesce  into  a  complex 
conjugate  pair  whose  real  part  crosses  zero  for  a  certain  value  of  z^.  Within 
the  accuracy  of  Figures  1  through  3,  the  dgenvalues  A  are  identical  to  those 
computed  fay  either  the  coupled  or  the  uncoupled  system,  or  the  perturbation 
equations  (18).  There  is  a  very  small  difference  for  the  surge  eigenvalue  as 


9 


shown  in  Figure  4,  but  again  the  agreement  between  coupled  and  uncoupled 
systems  is  very  good. 


C.  DEGREE  OF  STABILITY 

The  eigenvalues  of  the  reduced  order  system  (18)  designate  stability  or 
instability  of  motion.  A  single  measure  of  stability,  the  “degree  of  stability”, 
e«,  can  be  defined  as  the  maximum  real  part  of  all  three  eigenvalues  of  (18). 
This  measures  the  slowest  exponential  convergence  to  the  equilibrium  when 
negative  and  the  fastest  exponential  divergence  from  the  equilibrium  when 
positive.  Fbr  all  numerical  calculations  in  this  work,  the  degree  of  stability 
corresponds  to  the  real  part  of  a  complex  conjugate  pair  of  eigenvalues.  Typ¬ 
ical  results  are  presented  in  Figures  5  through  7.  Figure  5  shows  the  degree 
of  stability  versus  xgb  constant  dimensionless  speed  uq  =  0.5  and  different 
values  of  zqb'  Similar  results  are  shown  in  Figure  6  for  constant  xqb  =  0.015 
and  different  values  of  uo.  Finally,  a  three  dimensional  representation  is 
depicted  in  Figure  7. 


D.  CRITICAL  SPEED 

The  parameter  value  where  the  real  part  of  the  complex  conjugate  pair  of 
eigenvalues  shown  in  the  previous  figures  crosses  zero  defines  the  point  where 
linear  stability  is  lost.  This  critical  point  can  be  computed  by  considering 
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equation  (19).  Routh’s  criterion  ^plied  to  this  cubic  yields  A^Di  =  B^Ci 
which  can  be  solved  for  the  dimensionless  weight, 


where, 


W  = 


Ca.o  =  Zy,{M^  -  mxa)  -  M^{Z^  +  m)  , 
Ca.i  =  (m  -  Z,i,){zGB  cos  -  *ob  rin  ^o)  » 
Dj.i  =  Zit,{xoB  6fj  —  zqb  cos  flo) . 


(20) 


The  value  of  the  critical  speed  Ue  can  then  be  evaluated  from  (20)  end  (14). 
Typical  results  are  presented  in  Flgtues  8  and  9,  nondimensionalized  with 
respect  to  nominal  vehicle  speed  2.44  m/sec  and  length  4.26  m.  Vertical  plane 
motions  are  stable  for  forward  speeds  ki»  than  the  critical  speed.  It  can  be 
seen  that  stability  is  increasing  with  increasing  xg  while  xg  =  0  is  the  most 
conservative  condition  for  stability.  Therefore,  a  vehicle  which  is  stable  when 
properly  trimmed  will  remain  stable  for  off-trim  conditions.  Fbr  comparison, 
we  note  that  the  simple  stability  coefficient  defined  in  equation  (1),  is 
monotonically  decreasing  and  becomes  more  negative  for  decreasing  x^,  as 
shown  in  figure  10.  Thus  it  would  have  predicted  unstable  motions  for  the 
entire  range  of  parameters  shown  in  Figures  8  and  9.  For  completeness,  the 
value  of  the  steady  state  pitch  angle,  0o  (in  degrees),  is  shown  in  Figure  11 
versus  x^b  using  zqb  as  the  parameter.  This  pitch  angle  is  computed  from 
equation  (8). 
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Figure  11:  Steady  state  pitch  aagle  8q  Terras  mcB  for  zoB  fcotn.  0.00S  to  0.025  with 
increments  of  0.006 


m.  BIFURCATION  ANALYSIS 

A.  INTRODUCTION 

L/  all  cases  of  stability  loss  of  tbe  previous  chapter,  one  pair  of  com¬ 
plex  conjugate  eigenvalues  of  the  corresponding  eigenvalue  problem  crosses 
transversally  the  imaginary  axis.  A  situation  Hke  this  in  which  a  certain 
parameter  is  varied  such  that  the  real  part  of  one  pair  of  complex  conjugate 
eigenvalues  of  the  linearized  system  matrix  crosses  zero,  results  in  the  qrs- 
tem  leaving  its  steady  state  in  an  oscillatory  manner.  This  loss  stability 
is  called  Hopf  bifurcation  and  generically  occurs  in  either  supercritical  or 
subcritical  form.  In  the  supercritical  case,  stable  limit  cycles  are  generated 
after  the  nominal  straight  line  motion  loses  its  stability.  The  amplitudes  of 
these  limit  (^cles  are  continuously  increasing  as  the  parameter  distance  from 
its  critical  value  is  increased.  For  small  values  of  this  criticality  distance  the 
resulting  limit  <7cle  is  of  small  amplitude  and  differs  little  from  the  initial 
nominal  state.  In  the  subcritical  case,  however,  stable  limit  cycles  are  gener¬ 
ated  before  the  nominal  state  loses  its  stability.  Therefore,  depending  on  the 
initial  conditions  it  is  posrible  to  (fiverge  away  from  the  nominal  straight  fine 
path  and  converge  towards  a  limit  cycle  even  before  the  nominal  motion  loses 
its  stability.  This  means  that  in  the  subcritical  Hopf  bifurcation  case  the  do¬ 
main  of  attraction  of  the  nominal  state  is  decreasing  and  in  &ct  it  shrinks 
to  zero  as  the  critical  point  is  i^proached.  Random  external  disturbances  of 
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sufficient  magnitude  can  throw  the  vehicle  df  to  an  oscillatory  steady  state 
even  though  the  nominal  state  may  still  remain  stable.  After  the  nominal 
state  becomes  unstable,  a  discontinuous  increase  in  the  magnitude  of  motions 
is  observed  as  there  exist  no  simple  stable  nearby  attractors  for  the  vehicle 
to  converge  to.  Distinction  between  these  two  qualitatively  different  types  of 
bifurcation  is,  therefore,  essential  in  design.  The  computational  procedure 
requires  higher  order  approximations  in  the  equations  ci  motion  and  is  the 
subject  of  this  chapter. 


B.  THIRD  ORDER  EXPANSIONS 

The  nonlinear  heave/pitch  equations  of  nootion  (2),  (3),  and  (4)  are  writ¬ 
ten  in  the  form, 

9  =  q,  (21) 

w  =  Onto  +  oijq  +  ai3(«G!n  cos  0  +  xqb  an  +  d^(To, q) 

+ci(w,9)  ,  (22) 

q  =  Ojito  +  0339  +  «33(*GB  COS  +  xob  Mn  )  +  dl,(to,  q) 

+03(10,9)  ,  (23) 

where 

D„  =  (m-  -  Mi)  -  (mxo  +  Z^)(mxo  +  M^)  , 

oiiDv  =  (Jy  —  Mi)Z^  +  (mzG  +  Zi)M^  , 

013D*  =  {ly  -  Mi){m  +  Z,)  +  {mza  +  Z^)(Af,  -  m®c) , 
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<h*jD,  =  -(nuBo  +  Z^)W  , 

o>i\D^  =  (m  —  Z^)M^  +  (jtiatQ  +  , 

onD,  =  (m  —  —  mma)  +  (mco  +  M^){m  +  Z,)  , 

omD.  =  -(m-Z^)W, 

«i.(w,9)D.  =  (/»  -  +  (m*G  +  5^^)/,  , 

4,(w,9)Z),  =  (m-Z*)/,+(in*o+Af*)r^  , 

Ci(w»9)-0,  =  (In  —  M^)mxa^  —  {mma  +  Z^ymz^wq  y 
Ci(wyq)D,  =  -(m  -  Z^)mxav>q  +  (mco  +  M^)mxaq^  , 
and  7^,  7,  ate  the  cross  flow  integrals 

^  A(e)(to~«9)lti>~«9|de  ,  (24) 

^  idB  -  *«1*  ^  •  (25) 

The  system  of  equations  (21)  through  (23)  is  written  in  the  compact  form 

x  =  Ax+g(x),  (26) 

where 

X=[fl,ii;,g],  (27) 

is  the  three  state  tables  vector,  and  A  is  the  linearised  sytem  matrix  eval¬ 
uated  at  the  nominal  point  x*.  The  term  g(x)  contains  aU  nonlinear  terms 
ci  the  equations.  Hopf  Infurcation  analysis  can  be  performed  by  isiflating 
the  primary  nonhnear  terms  in  g(x).  Keeping  terms  up  to  third  order,  we 
can  write 

g(x)  =  gt»)(x)+gt»)(x).  (28) 
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Using  equations  (21)  through  (25),  the  various  terms  in  (28)  can  be  written 
«•> 

-  n 

9\  —  0 , 

=  (/y  -  Mi)mxQt^  -  {mxa  +  Z^ymxotoq  +  ^^{v},q)  ,  (29) 

=  -(m  -  Z^)m*a^q  +  (m*©  +  +  <i^\v>,q)  , 

and 

9?^  =  0, 

=  d^^te,  q)  +  Jois(*oB  rin  ^  -  sqb  cos  tfo)^  ,  (30) 

^*)  -  <^*)(ty ,  q)  +  Jaas(xoB  rin  8s  -  ^B  cos  8i>)8*  . 

Expansion  in  Thylor  series  of  d^,  4,  requires  expansion  of  the  cross  flow 
integrals  /«,/,,  which  require  the  IkylOT  series  of 

m = ^1^1 .  (31) 

This  expression  can  be  converted  into  an  analytic  function  using  Dalxell’s 
approximation  (Dalzell,  1978), 

which  b  derived  by  a  least  squares  fit  of  an  odd  series  over  some  assumed 
range  of  namely  — <  (  <  (e*  appnndmation,  which  b  shown  in 
figure  12,  has  been  extensively  used  in  ^p  roll  motion  studies  and  b  very 
useful  for  its  intended  purpose.  However,  in  the  present  problem  it  suffers 
from  several  major  drawbacks: 


26 


•  It  introduces  a  finear  term  which  depends  on  the  assumed  range  of 
motion,  and  it  renders  the  critical  speed  hmction  ai  the  vehicle  motions. 

a  The  cubic  term,  which  is  ultimately  responsible  for  the  Hopf  bifurcation 
analysis,  is  a  function  of  the  assumed  range  of  vehicle  motions  which 
can  not  be  known  in  advance. 

•  As  Figure  12  demonstrates,  the  dope  of  the  actual  curve  at  the  origin 
is  dgnificanly  different  than  the  approximation,  which  would  make  the 
bifurcation  results  unreliable. 

Instead  of  Dalzell’s  approximation,  we  employ  the  concept  of  generalized 
gradient  (Clarke,  1983),  which  is  used  in  the  study  of  control  qrstems  in* 
volving  discontinuous  or  non-smooth  functions.  In  this  way  we  approximate 
the  gradient  of  a  non-smooth  fimction  at  a  discontinuity  a  map  equal  to 
the  convex  dosure  of  the  limiting  gradients  near  the  discontinuity.  In  our 
problem  we  write, 

/({)  =  &I6I  +  2|<ol«  -  6.)  +  -u?  +  f'HO  .  (33) 

as  the  Taylor  series  epansion  of  f(()  near  (p*  The  dgn  function  in  (33)  can 
be  approximated  by, 

dgn((o)  =  l^tanh  .  (34) 

A  graphical  representation  of  the  approximation  (34)  is  shown  in  Figure 
13.  The  quantity  7  is  a  small  regularization  parameter  and  is  tised  in  the 
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Figure  12:  Graphical  repr Mentation  of  Dalseil*s  apprarimation  of  |  vartns  f/fe- 
Solid  corre  is  the  exact  expression  end  dotted  carve  is  the  apprwnritnat^nn  (33) 
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next  section  for  the  proper  normalisation  ct  the  results.  Using  (34),  we  can 
approximate  /(()  in  the  vicinity  of  ^  =  0  by, 

<l<l  «  •  (3») 

Since 

(>-^w-zq  ,  (36) 

we  can  express  the  non-smooth  cross  flow  integral  terms  by, 

J.  =  ^(E,yi‘-3B,i.i‘q  +  3E,w^-E,f),  (37) 

Iq  —  —  £1^)  ,  (38) 

67 

where 

Ei  =  z*b{m)  dm  ,  (39) 

are  the  moments  of  the  vehicle  ‘Sraterplane”  area. 

C.  COORDINATE  TRANSFORMATIONS 

Using  the  previous  second  and  third  OTder  Taylor  series  expanrions,  equa¬ 
tion  (26)  is  written  in  the  form, 

X  =  Ax  -I-  •  (40) 

If  T  is  the  matrix  of  eigenvectors  of  A  evaluated  at  the  critical  point  u  = 
the  finear  change  of  coordinates, 

X  =  Ts  ,  B  =  T'‘x  ,  (41) 
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Figure  13:  (baphical  lepreeenttion  of  the  rign  fimctkm  end  iti  hjrperbolie  tengcBt 
ai»prcaimstioii  (34).  S(did  corre  oocrespeads  to  7  s  0.1  and  dotted  eorre  oocxe- 
^onds  toy  —  0.01 


transfonns  qrstem  (40)  into  its  normal  oooniinate  form, 

a  =  T"*  ATa  +  T"V*^(Ta)  +  T-V*)(Ta)  .  (42) 

At  the  Hopf  bifurcation  point,  matrix  T~^AT  takes  the  form, 


T-^AT 


0  -a»o  0 

(Ub  0  0 

0  0  p 


where  uo  is  the  imaginary  part  of  the  critical  pair  of  eigenvalues,  and  the 


remaining  eigenvalue  p  is  negative.  Fbr  values  of  u  dose  to  the  bifurcation 


poit  He,  matrix  T~^AT  becomes, 


T-'AT 


o/e  — (a»o  +  u/e)  0 

(ufo  +  o/e)  o^€  0 

0  0  P  +  P^e 


where  e  denotes  the  criticality  difference 


e  =  U  —  Ue  , 


(43) 


and 

a'  =  derivative  of  the  real  part  of  the  critical  eigenvalue 
with  respect  to  e  , 

u/  =  derivative  of  the  imaginary  part  of  the  critical  eigenvalue 
with  respect  to  e  , 

jf  =  derivative  of  p  with  respect  to  e . 
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D.  CENTER  MANIFOLD  EXPANSIONS 

Due  to  continuity,  the  eigevalue  p+j/e  remains  negative  for  small  imncero 
values  of  e.  Therefore,  the  coordinate  zg  corresponds  to  a  negative  eigenvalue 
and  is  asymptotically  stable.  Center  manifold  theory  predicts  that  the  rela¬ 
tionship  between  the  critical  coordinates  ,  S3  and  the  stable  coordinate  Zg 
is  at  least  of  quadratic  order.  We  can  then  write  zg  as, 

Zi  =  aiis^  -h  oitgZiZj  -I-  a33s|  ,  (44) 

where  the  coefRcients,  Oij,  in  the  quadratic  center  manifold  expansion  (44) 
need  to  be  determined.  By  differentiating  equation  (44)  we  obtain, 

za  =  2aiizizi+aig(ziz2+ziz2)-h2a„z2za  .  (45) 

We  substitute  ii  =  — W0S3  and  zg  =  u>o*i  from  equation  (42)  into  (45),  and 
we  obtain 

is  =  aiaWto^  +  2(0,3  -  an)wo*i*a  -  •  (46) 

The  third  equation  of  (42)  is  written  as, 

i,=p»,+  [T-*g<’>(Tm)]|^l  ,  (47) 

where  terms  up  to  second  order  have  been  kept.  If  we  denote  the  dements 
of  T  and  T“^  by, 

T  =  [rn,,]  ,  T-‘  =  K)  ,  (48) 


then 


fdx 

T-V’HT*)  =  dg 


7 


.jMmmm}:  iipiii 


where 

di  =  ni2(/n^  +  +  ^a7*i)  +  +  ^*i*a  +  ^7*j) ,  (49) 

dj  =  +  ^36^1^  +  ^»7*j)  +  TJa^C^a#*?  +  ^»*i*7  +  far^)  »  (50) 

4s  =  7*3a(/a«*i  +  ^26*1^  +  ^27*j)  +  7»3s(^M^  +  ^88*1  *J  +  ^7^)  •  (51) 

Expressions  for  the  coefficients  iij  are  (pven  in  the  Fbrtran  prc^ams  in  the 
^pendix. 

Equation  (47)  then  becomes 

is  =  pzs  +  dt  ,  (52) 

and  substituting  (44)  and  (51)  into  (52)  we  get, 

is  =  (pan  +  Tisa^s  +  7»33^)X|  +  (pan  +  tisj/m  +  713s4m)^i^ 

+(pa3a  +  nsisln  +  7i9s^7)xa  .  (53) 

Comparing  coefficients  of  (46)  and  (53)  we  get, 

-pan  +  <*»oaia  =  tjsa^ii  +  tiss^sb  ,  (54) 

— 2ci;oaii  —  paia  +  2u;oaaa  =  nsiln  +  tim/m  »  (55) 

— o/Qaia  —  poaa  =  fhaitr  +  tiss^  •  (56) 

Solution  of  the  system  of  linear  equations  (54)  through  (56)  yields  the  coef¬ 
ficients  in  the  center  manifold  expansion  (44). 
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E.  AVERAGmC 


Using  the  previous  Tsylor  expansions  and  center  manifold  approxima¬ 
tions,  we  can  write  the  reduced  two-dimensional  system  that  describes  the 
center  manifold  flow  of  (42)  in  the  form, 

ii  =  a'exx  -  (a/o  +  w'«)*a  +  Fi(xi,arj)  ,  (57) 

ij  =  (wo  +  o>*€)zi  -h  a'e*,  +  F2(xi,zj)  ,  (58) 

where 

■fi(*ii*a)  =  ♦*11^  +**ia*i*i  +Titziz^  +**14^^ 

+fhi*i  +Pia*i*a  +Pi8*j  »  (59) 

^a(*ii*a)  =  T’ai^  +  i*aa^*a  +  »*as  •*!*»+ »*a4*s 

+l>aj*?  +  Pa8*i*a  +  Pn4  »  (®0) 

and 

ry  =  riia/jj  -f-  i  t  =  1,2  ,  j  =  1, . . .  ,4  ,  (61) 

Pij  =  ^Ka^afc  +  n<s^  »  *  =  1, 2  ,  j  =:  1, 2, 3  ,  ks=  j  +  4  .  (62) 

we  introduce  polar  coordinates  in  the  form, 

zi=iZcos^,  zj^Asin^,  (63) 

we  can  use  (57)  and  (58)  to  produce  an  equation  describing  the  rate  of  change 
of  the  radial  coordinate  R, 


R  =  a’eR  +  +  Q{^)R^  • 


(64) 


This  equation  contains  one  variable,  12,  which  is  slowly  varying  in  time,  and 
another  variable,  which  is  a  fast  variable.  Therefore,  equation  (64)  can  be 
averaged  over  one  complete  cycle  in  ^  to  produce  an  equation  with  constant 
coefficients  and  similar  stability  properties, 

R  =  a'eR  +  KI^  +  LR*  ,  (65) 

where 

K  =  ^  Jq  ^  ~  »  (6®) 

L  =  =  (67) 

Therefore,  the  averaged  equation  (65)  becomes 

R  =  a'tR  +  KH*  .  (68) 


F.  LIMIT  CYCLE  ANALYSIS 

Equation  (68)  admits  two  steady  state  solutions,  one  at  i2  =  0  which 
corresponds  to  the  trivial  equilibrium  solution  at  zero,  and  one  at 

i2o  =  •  (69) 

This  equilibrium  solution  corresponds  to  a  periodic  solution  or  limit  <^cle  in 
the  cartesian  coordinates  zi,  Z3.  For  this  limit  cycle  to  exist,  the  quantity  i2o 
must  be  a  real  number.  In  our  case  ci  is  always  positive,  since  the  system 
loses  its  stability;  i.e.,  the  real  part  of  the  critical  pur  of  rigenvalues  changes 
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from  negative  to  pontive,  for  increasing  u.  Therefore,  existence  of  these 
periodic  solutions  depends  on  the  value  of  K.  Specifically, 

•  if  iiT  <  0,  periodic  solutions  exist  for  e  >  0  or  u  >  ««,  and 
s  if  IC  >  0,  periodic  solutions  exist  for  e  <  0  or  it  <  Ue> 

The  characteristic  root  of  (68)  in  the  vicinity  of  (69)  is 

)9=-2a'e,  (70) 

and  we  can  see  that 


•  if  periodic  solutions  exist  for  u  >  Ue  they  are  stable,  and 

•  if  periodic  solutions  exist  for  it  <  Ue  they  are  unstable. 


The  period  of  these  periodic  solutions  can  be  estimated  as  follows.  Equa¬ 
tions  (57),  (58),  and  (63)  produce  an  equation  in  ^  similar  to  (64), 

^  =  wo+u,'€  +  -H  Gi<l>)R  .  (71) 

The  averaged  form  of  (71)  is 

i  =  wo+uj'€  +  MR^  ,  (72) 


where 

=  (73) 

The  limit  cycle  period  can  be  computed  by  substituting  (69)  into  (72), 


Mo  -I-  Mi2o 
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G.  RESULTS  AND  DISCUSSION 


Typical  results  ci  the  nonlinear  stability  coefficient  K  are  shown  in  Fig¬ 
ures  14  and  15.  Figure  14  presents  a  {dot  of  K  -y  versus  cq  for  so  =  0.015 
and  for  different  values  of  Cj).  It  shotild  be  emphasised  that  the  use  ci  K  -y 
is  more  meaningful  than  the  use  of  FC,  mnce  it  {voperly  accounts  for  the 
use  of  the  regularization  parameter  y  as  seen  from  equations  (35)  and  (68). 
Numerical  evidence  demonstrates  that  all  curves  K  •  y  versus  ma  converge 
for  y  —*  0.  For  practical  purposes,  values  of  y  smaller  than  0.001  produce 
identical  restilts.  The  results  of  Figure  14  demonstrate  the  profound  effect 
that  the  quadratic  drag  coefficient  Cd  has  on  stability  of  Emit  cycles.  All 
Hopf  bifurcations  are  supercritical  (K  <  0),  and  they  become  stronger  su¬ 
percritical  as  C7o  is  increased.  It  is  worth  noting  that  results  for  Cd  =  0 
produce  subcritical  behavior,  FT  >  0,  which  is  clearly  incorrect.  Thus,  ne- 
(^ecting  the  effects  of  Cd  would  have  produce  entirely  wrong  results  in  the 
present  problem.  Figure  15  shows  a  plot  of  FT  •  7  versus  xq  for  Cd  =  0.5 
and  different  values  of  the  metacentric  height  zc.  It  can  be  seen  that,  the 
bifurcations  become  stronger  supercritical  as  initial  stabiEty  xq  is  increased. 

The  bifurcation  analysis  results  are  verified  by  direct  numerical  simula¬ 
tions  shown  in  Figures  16  through  18.  Figure  16  shows  the  results  of  two 
numerical  simulations  for  two  values  of  nominal  speed  vo  in  terms  of  the 
vehicle  pitch  angle  B  (in  degrees)  versus  time  (in  seconds).  The  critical  value 
of  speed,  tie,  is  about  0.495  as  can  be  seen  from  Figure  8,  while  the  other 
parameters  in  the  simulations  were  Cd  =  0.5,  »g  =  0.015,  and  xa  =  0.  It 


37 


1^01 _ I _ I _ I  ‘ _ I _ I - 1 - 1 - 1 - 

*_10  -e  -6  >*4  -2  0  2  4  6  B  10 

xl0“» 


Figure  14:  K  -y  versus  Xfj  iat  zq  =  0.015  and  different  ndues  of  the  drag  codhcient 

Co 
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FIgaxe  15:  £*  •  7  vertot  zq  fm  Cd  —  0.5  and  different  valnes  of  the  metacentric 
height  XQ 
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Figure  16:  Tinie  hiitoiies  for  Co  =0^,Jta  =  0.015,  cg  »  0,  and  two  difirat 
values  of  nomiiial  speed 
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am  be  teen  that  conTexgence  to  aero  ii  entured  fi>r  tio  <  «■«  and  a>nvergence 
to  a  fimit  cycle  occurs  fi>r  110  >  ««•  This  indicates  supercritical  behaTior  as 
shown  before.  A  selection  of  time  histories  is  shown  in  Figure  17  for  a  range 
of  forward  speeds  and  the  same  parameters  as  in  Figure  16.  The  same  initial 
disturbance,  0  =  5  degrees,  was  introduced  at  <  =  0  for  all  simulations.  It 
can  be  seen  that  the  amplitude  oi  limit  cycles  increases  as  the  distance  of  it 
from  lie  is  increasing.  The  rate  of  convergence  of  solutions  to  their  limit  cy¬ 
cles  is  also  increasing,  while  their  period  remains  essentially  constant.  These 
results  are  sum  irised  in  Figure  18,  where  the  amplitudes  of  the  numeri¬ 
cally  computed  limit  cycles  are  (dotted  versus  uq.  The  behavior  is  clearly 
supercritical,  which  agrees  with  our  findings  of  the  bifurcation  analysis. 


43 


IV.  BIAS  EPPECTS 


A.  LOSS  OP  STABILITY 

Stability  analysis  ci  motions  at  a  nonsero  angle  o£  attack  can  be  per¬ 
formed  by  first  introducing  some  bias  into  the  steady  state  solution  and  its 
perturbations.  This  can  be  achieved  by  maintaining  a  nonsero  dive  plane 
angle  at  nominal.  In  this  case  the  steady  state  solutions  are  ^  =  0,  and  to^, 
6o  are  computed  from, 


Z<^voo  ~  Cn -Eotoo  |too  I  Z$6  =  0  ,  (75) 

M^wq  +  Co-Sriti)io|ti>o|  ~  oos^  +  +  ACjtf  =  0  .  (76) 


The  coefficients  JSq,  Ei  are  computed  from  (39).  hi  order  to  solve  (75)  we 
observe  that  when  ^  >  0,  then  to^  <  0, 


Wo  = 


-Z^-y/Zl  -ACoEoZfS 

-2CdEo 


(T7) 


and  when  S  <0,  then  too  >  0, 


— Zw  —  yjz^  +  4CdEoZsS 

-2Cofi> 


(78) 


The  angle  6o  can  then  be  computed  from  equation  (76). 


Linearisation  of  the  equations  of  motion  in  the  neighborhood  of  the  above 
equilibrium  point  produces  the  linear  system, 


(t7»  —  Z^yd}  —  {rnzQ  +  Z^)9  —  (Z*  —  2CuEo^ol}^ 
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+(Z,  +  m  +  2Co£ri|iiHi|)4  t 


(79) 


+  msBa)w  +  (/,  —  M^)i  —  (Af*  +  IwoDw 
+(Jli,  -  m*o  -  nuoiflo  —  2Ci>£9|ti>o|)9 
+fr(«oB  m  tfij  -  zoB  CO*  ,  (80) 

8  =  g,  (81) 

where  the  variables  w,  8  are  understood  as  small  deviations  from  their  equi¬ 
librium  values.  Numerical  solution  of  the  generalised  eigenvalue  problem  (79) 
through  (81)  yields  the  critical  speed  values  where  the  nominal  equilibrium 
solution  becomes  unstable. 

B.  ANALYSIS  OF  HOPF  BIFURCATIONS 

It  can  be  numerically  verified  that  the  above  calculations  for  the  new 
critical  speed  result  in  a  loss  of  stability  in  tiie  form  of  Hopf  bifurcations,  as 
for  the  8  =  0  case.  These  Hopf  bifurcations  can  be  analysed  using  the  same 
general  methodology  that  was  developed  in  the  previous  ch^ter.  The  non¬ 
linear  expansions  are  like  equations  (21)  to  (^)  vrith  the  follovring  changes: 
The  substitutions 

^  —  2Cx7£otU|)  , 

-f-  2Cn^ioo  , 

Z^  *}■  rn  — ►  Z^  tn  -H  2C!oEiV>q  , 

M,  —  msDa  —*■  M,  —  mxo  —  mzawo  —  2CDE2too  , 
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are  assumed  in  the  definition  of  coefficients  Oij.  FVirthemmre,  equation  (33) 
prduces  2nd  order  contributions  due  to  wq  0,  as  well  as  3rd  order.  Using 
(33)  we  can  compute  the  2nd  order  expansions  of  and  in  (29)  using, 

i<?>  =  CD(EoV}^-2Erwq  +  Etq^),  (82) 

=  Cd{,E\w^  —  2E^wq  +  Bl^^) .  (83) 

These  equations  are  valid  for  wq  >  0  or  f  <  0.  For  £  >  0  the  signs  of  Ek 
must  be  switched.  The  third  order  expansions  and  are  the  same 

as  in  equations  (37)  and  (38).  Using  these  additional  terms,  the  nonlinear 
stability  coefficient  K  can  be  computed  in  the  same  way  as  in  the  previous 
chapter. 


C.  RESULTS  AND  DISCUSSION 

Typical  results  for  ^  ^  0  are  shown  in  Figures  19  through  24.  Figure 
19  shows  the  equilibrium  pitch  angle  ^  (in  degrees)  versus  zq  for  different 
values  of  S  from  —5  to  5  degrees  with  increments  of  1  degree.  Solid  curves 
correspond  to  positive  S  and  dashed  curves  to  negative.  Figure  20  shows  the 
degree  of  stability  for  the  equilibrium  points  of  Figure  19,  while  Figure  21 
presents  the  degree  of  stability  in  a  three  dimensional  view.  It  can  be  seen 
that  positive  and  negative  values  of  6  have  almost  identical  stability  charsu:- 
teristics.  Furthermore,  the  degree  of  stability  becomes  more  negative  as  the 
absolute  value  of  £  is  increased,  which  means  that  we  expect  a  wider  domain 
of  stability  in  this  case.  This  is  verified  by  the  critical  speed  {dots  shown  in 
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Figure  19:  Pitch  an|^e  ^  vemu  so  fiar  s  0.015,  uo  —  0.5,  and  d 
atS 
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figure  20:  Degree  of  stability  tosus  zq  ^  xq  ^  0.015,  uo  s  0.5,  and  <fiffereiit 
’values  of  B 
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degree  ot  stability 


Flgoie  21:  Degree  of  stability  verstu  xq  and  uq  fi>r  s^  =  0.015  and  two  ^adaes  of  S 
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Figoie  22:  Critical  spaed  tie  vamu  zq  bt  zo  ~  0.015  and  different  ndues  of  S 
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Figure  24:  Nonlinear  stabilitj  coefficient  K  versus  xg  fat  Co  =  0.5,  zq  =  0.015, 
and  different  values  of  6 
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Figures  22  and  23.  The  critical  speed  is  minimuin  for  S  ^0  and  it  increases 
monotonically  with  increasing  absolute  value  of  8.  This  stabilising  effect  of 
asymmetry  (bias)  remains  ^proximately  true  for  the  nonlinear  analysis,  as 
demonstrated  the  results  of  Figure  24.  It  can  be  seen  that  the  nonlinear 
stability  coefficient  K  becomes  more  negative  as  ^  is  decreasing  From  sero. 
For  increasing  8,  K  becomes  less  negative  but  the  difference  From  the  £  =  0 
calculations  appears  to  be  very  small.  Therefore,  limit  cycle  stability  is  not 
agnificantly  affected  by  the  bias  effects  that  are  induced  by  small  nonsero 
dive  plane  angles. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  thesis  presented  a  comprehensive  nonlinear  study  of  straight  line  stabil¬ 
ity  of  motion  of  submersibles  in  the  dive  plane  under  open  loop  conditions. 
A  systematic  perturbation  analysis  demonstrated  that  the  effects  of  surge  in 
heave/pitch  are  small  and  can  be  neglected.  Primary  loss  oi  stability  was 
shown  to  occur  in  the  form  of  Hopf  Infurcations  to  periodic  solutions.  The 
critical  speed  were  instability  occurs  was  computed  in  terms  metacentric 
height,  longitudinal  separation  of  the  centers  of  buoyancy  and  gravity,  and 
the  dive  plane  angle.  Analysis  of  the  periodic  solutions  that  resulted  from  the 
Hopf  bifurcations  was  accomplished  through  Taylor  expansions,  up  to  third 
oxdcTf  of  the  equations  of  motion.  A  consistent  approximation,  utilising  the 
generalized  gradient,  was  used  to  study  the  non-analytic  quadratic  cross  flow 
integral  drag  terms.  The  results  indicated  that  loss  of  stability  occurs  always 
in  the  form  of  supercritical  Hopf  Infurcations  with  stable  limit  cycles.  It  was 
shown  that  this  is  mainly  due  to  the  stabilizing  effect  of  the  drag  forces  at 
high  angles  of  attack.  As  a  recommendation  for  further  research,  we  suggest 
a  nonlinear  analysis  of  coupled  motions  in  aiz  degrees  of  freedom. 
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APPENDIX 


The  following  is  a  lilt  and  description  of  the  computer  programs  used  in  this 
thesis.  The  programs  are  written  in  FORTRAN  or  MATLAB.  Complete 
printouts  of  the  programs  follow  after  the  list. 

a  PERTURB.M 

MATLAB  program  for  performing  surge/heave/pitch  perturbation  anal¬ 
ysis. 

a  CR1T-0.M 

MATLAB  program  for  calculating  the  critical  speed  for  j  =  0. 
a  GRIT  J}ELTA.M 

MATLAB  program  for  calculating  the  critical  speed  for  £  ^  0. 
a  HOPF.O.FOR 

FORTRAN  program  for  calculating  the  nonlinear  stability  coefficient 
K  for  6  ^0. 

a  HOPF  J>ELTA.FOR 

FORTRAN  program  for  calculating  the  nonlinear  stability  coefficient 
KhtS^O. 

a  SIM.FOR 

FORTRAN  program  for  simulating  the  equations  of  motion. 
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t  Pregraa  partnrb.a 

X  Util  fil«  iqppllM  eoBe«p«a  of  portnrbotion  analysia  thoory, 
X  in  ordor  to  provo  that  tko  3  by  3  lystMi  doacribing  aotioa 
X  in  tho  Tortienl  piano,  doeonploa  frea  tbo  4  by  4  ono. 


zbo  ■ 

1.94; 

6  ■ 

32.2: 

L  ■ 

13.9792; 

Bdi  - 

0.6orho*L*2: 

lid2  ■ 

0.6orho*L*3: 

iid3  - 

0.6*rho*L'*4: 

nd4  ■ 

C.BarhooL'S; 

iy  - 

661.32/nd4: 

aqdot 

■  -6.330-4; 

sadot 

-  -1.46290-2; 

*q 

-  7.6460-3; 

sa 

•  -1.3910-2; 

■qdot 

■  -8.80-4; 

■adot 

■  -6.610-4; 

■  -3 .7020-3; 

m 

■  1.03240-2; 

od 

■  0.016; 

A 

•  0/L; 

■ 

■  1666.2363/ 

zndot 

■  -0.06oa; 

xb 

-  0/L; 

no 

■  4; 

a 

-  1666.2363. 

b 

■  a; 

>8 

•  0; 

>8 

■  linapaco(0 

xgb  ■ 

xg-xb; 

agb  ■ 

xg-ab; 

thota 

■  ataa(-xgb./ 

alpha 

■  ni-zadot: 

bota 

■  -za; 

Ai  ■  a-xndot; 

B1  - 

-2oed; 

.”2); 
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A  ■  Oi-svde«)«(ij-»qdet)-(B«xg«aqdot)*(Mide««s*zg); 

-(B«xg«sqdot ) *am ; 
for  i  ■  l:longth(sg) 

doltu(i)  ■  xgb*«««ia(tliota(i))-xgb(i)*v*ee«(thotm(i)); 

C(i)  ■  -dolt««(i)*(B-**dot)+aw*(Bq-B**g)-(*q‘**)*wi; 

0(1)  ■  doltas(i)«x«; 

A2(l)  ■  (»**g<i))'2*(B-«»dot); 

B2(l)  -  -(■•«g(i))“2*«w; 

CA  -  >A1*A: 

C8  ■  -A1*B'*'B1«A: 

CC  -  -Al*C(i)^Bl*B: 

CD  -  -Al*D(i)+Bl*C(i) ; 

CS  -  Bi*0(i) : 

p  -  CCA  C8  OC  CD  CE]; 
r(l:4.i)  -  zoots(p); 

l«aida(l:4,i)  ■  r(l;4,i)-(r(l;4,i)  .*3.*(»lpl«i.*r(l;4,i)+bot%)). 

.♦(■.•*g(i)) .*2./(4.*CA.*r(l:4,l) .“3+. . . 
3.*CB.*r(l:4,i).*2+2.*CC.*r(l:4,i)+CD); 

PA  -  -A1*A+A2(i}: 

PB  -  -AloB-^BloA+BBCi) ; 

PC  •  -Al*C(i)+Bl*B; 

PD  -  -Al*0<i)+Bl*C(i); 

PE  -  B1*D(1): 

p4  ■  CPA  PB  PC  PD  PE];  X  o-raluoa  of  4  by  4  aystoiB. 
x4(l:4,i)  ■  root>(p4); 
ond 
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t  Progm  erlt.O.B 

t  Evaluation  ot  critical  ^ood  for  dalta"0 

rho  ■  1.94; 
g  >  32.2; 

L  -  13.9792; 
ndl  ■  0.6*t1io*L*‘2; 
nd2  ■  0.6*rlio«L"3; 
nd3  *  0.6«xlio*L'‘4; 
nd4  ■  0.6*rhe*L'‘6; 
i7  -  661.32/nd4; 
zqdot  <■  -4.33a*4; 
zvdot  -  -1.4629«-2; 
zq  ■  7.64Sa~3; 

ar  ■  -1.391«*2; 
nqdot  ■  -8.8o<-4; 
avdot  ■  -6.81a*4; 

■q  ■  -3.702a-3; 
m  •  1.0324O-2; 

cd  ■  0.016; 

zb  ■  0/L; 

m  •  1568. 2383/ (g*iid2); 

xndot  ■  -0.06am; 
xb  >  0/L; 

zg  >  linapaco(-0. 01, 0.01,40); 
zg  ■  liaspaco(0. 005, 0.026 ,40); 

Ov  ■  1  -  av.*(zqam)./(zw.a(mq-m.*zg)); 
lambda  ■  -2*cd/m-zndot ;  %  via  •  -1.8439 

zgb  ■  zg-zb; 
zgb  B  zg-zb; 

for  j-l:longtb(zg) 
for  iBl:longth(zg) 

thota(i,J)  ■  atan(-zgb(i)/zgb(j)); 

aO  ■  (m-zvdot)*(i7-mqdot)-(mvdot+m*zg(i))*(zqdotam*zg(i)); 
bO  ■  (-zndot«m-m*mv-zq*m)*zg(i)4(-m*mqazwdot*mq-zqdot4«r... 
-zq^Bvdot-m<*anrdot-i7*zwamqdot«zw) ; 
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■  >  JftiL  58^,  Jt 


cO  •  ~ai*sv*xg(i)-HBq<*n-aq«m-a*Mr: 

cl  ■  (-■*3tg(i)+»»4ot*xg(i)+m*xb-»«dot**b)*«ia(tk«ta(l,j)). . . 

♦(-■•xb-xwdot*xg(  j  )+xwdot*ab^*»g(  j»*eo.(th*ta(i ,  j  )  ) ; 
dl  -  (zw*xg(i)<-z«*xb)*«ia(tb«ta(i.j))«... 

(xw*xb-xw*xg(  j  )  )  *eo«(tliata(i ,  j  )  ) ; 


X  iltar  iqBplying  ib-BC... 


«(i,j)  •  bO*eO/(aO*dl-bO*el): 

tiO(i.j)  -  (l8S6.2363/(Bdl**(i,j)))*.5; 

cad 

aad 
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X  Fil«  erit .delta. ■ 

X  Evaluation  of  critical  apaad  for  aonsaro  atom  plana  angle. 

1.94; 

32.2; 

13.9792; 

0.6arho*L'‘2; 

0.6*rlio*L‘‘3: 

0.6arlio*L*4: 

0.6arho*L'5; 

Mi.32/nd4: 

■  -6.33a-4; 

-  -1.4S29a-2; 

-  7.546a-3: 

■  -1.391a*2: 

■  -8.8a-4; 

-  -6.6ia-4: 

-  -3.702a-3; 

>  1.0324a>2: 

•  1666. 2383/ (g*nd2); 

■  -0.06an; 

•  0/L; 

-  0/L; 

«  0.016; 


odz 

-  0.6; 

X  if 

edz"0,  division  by  zero  provides  Mall 

zd 

■  -6.603e-03; 

nd 

«  -2.4090-03; 

EO 

-  0.1036609; 

X  EO 

■■  Litegral  of  b(z}*dz 

El 

-  -2.36299820-03; 

X  El 

■>  Integral  of  z*b(z)*dz 

E2 

-  7.36721470-03; 

X  E2 

■■  Integral  of  z-2*b(z)*dz 

*8  ■ 

linzpaceC-O.Ol.O.Ol.OO) ; 

X  HOI-DIHEESIOIAL 

uO  ■  linspaea(2,6,40); 

xgb  •  xg-zb; 
zgb  •  zg'zb; 

vl  >  1566. 2383. /(ndl. *00.-2); 

delta  ■  input (’Enter  the  value  of  delta:  delta  ■  liaspaee(  ,  ,  )*); 
delta  ■  delta*pi/180; 


rho  ■ 

8  ■ 
L  ■ 

ndl  ■ 
nd2  - 
nd3  ■ 
nd4  ■ 
iy  - 
nqdot 
zvdot 
zq 
zv 

nqdot 

nndot 

nq 

a 

zudot 

zb 

zb 

*8 
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for  k>l:l«agtli(d«ltm) 
for  jal:loiigth(nO) 
for  i«l:loagtli(zg) 

if  dolta(k)  >  0 

«0  ■  -(s«-»oqrt(zv‘'2-4*eds*EO*ad*dolta(k)})/(2*eds*EO) ; 
eooffl  "  OMi*«0-eds*«0‘‘2«Bl‘Hkl*dolta(k))/«l(j); 
polyl  ■  [(xgb(i)“2+«gb*2)  (-2*eooffl*sgb)  (cooffl''2-rfb(i)'2)] 
thotaO  ■  oaiaCroots (polyl)): 

cboekl  ■  sgb(l)«eos(thotaO(l))-»^Sgb*si&(thotaO(l)); 
ehoek2  ■  zgb(i)«eot(thota0(2))-»agb*sia(thota0(2)) ; 
if  abs(ehock2-eooffl)  <  abaCekoekl-eooffl) 
tbotaO  »  tbota0(2): 

also 

thotaO  ■  tkotaO(l): 
oad 


B»C  Oa-awdot)  -(B*zg(i)’»xqdot)  0 
-(■sdot-Mi*zg(i))  iy-aqdot  0 
0  0  13; 

A"C(xw'f2aeds*E0aw0)  (xq<»m)-2*eda*Ei*wO  0 

Hi-2*edz*El««0  (■q-m«xg<i))-B*sg*«0424>eds*E2*«0  . . . 

(zgb(i) *aiii(thotaO) -zgb^eoa (tbotaO) )*vl  ( j ) 

0  1  0] ; 
olaoif  dalta(k)  <  0 

«0  >  -(zvfaqrt(sw''2-«-4*eda*E0*sd*dolta(k)))/(-2*edB4'E0); 
eooffl  ■  (jnr*wO+edaoirO“2*El+nd*dolta(k))/wl(j) ; 
polyl  ■  [(xgb(i)*‘2+*gb"2)  (-2oeoofflaxgb)  (coofll*‘2-xgb(i)*2)] 
thotaO  -  aainCrootaCpolyl)): 

ehoekl  ■  rgb(i)*eos(thotaO(l))<»^sgbasia(tbotaO(l)) ; 
eboek2  >  xgb(l)*eos(tliota0(2))'*'agb*sia(thota0(2)); 
if  abs(ehoek2'eooffl)  <  abs(eboekl*eo#ffl) 
thotaO  ■  thota0(2); 

olao 

thotaO  >  thotaO(l): 
and 

B»C  (m-aadot)  -(m*xg(i)+aqdot)  0 
-(B«dot<»n*xg(i))  iy-aqdot  0 
0  0  1]; 
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A«C(zv**2*eds*EO*wO)  (zq-nB)'*’2*ects«Sl*vO  0 

BW42*edz*El*«0  (Bq>a*zg(i))'m*zg««0-2*edz*E2«v0  ... 

(zgb(i)  •■iii(th«t«0)  *sgb*eos  (thataO)  )  *vl  (  j  ) 

0  1  0]; 

«nd 

avals  ■  aigd.B): 
dagstab(i.j)  ■  nazCzaalCavals)) : 

and 

and 

%  Evaluation  of  tha  critical  spaad.  ucrit: 

for  l>l:laagtb(xg) 
for  j«l:(laagtb(uO))-l 

flagl  ■  8ign(dagstab(i,j)): 
flag2  ■  8ign(dagstab(i,j-»l)) : 
if  flag!  '>  flag2 

uerit(i,k)  »  <uO(j+l)adagstab<i,j)-uO(j)*dagstab(i, j+1))/. . . 
(dagstabCi ,  j  ) -dagstabd ,  j  4 1 )  ) ; 
and 
and 
and 
and 

%  Ravalnation  of  tba  nominal  angla,  'thatacr' 

%  at  tba  critical  spaad,  ’ucrit’: 

%  Nota  that  in  this  casa,  at  a  eartain  ’xg’ , 
t  corraspends  a  spaeific  ’ucrit*. 

%  That '8  idiy  «a  usa  tha  sana  indax  *i*: 

«2  -  1666.2363./(iidl.*uerit.‘2): 

for  i"i:laagth(xg) 

if  dalta(l)  >  0 

«0  ■  -(zv4sqrt(zv*2-4*cdz*E0*sdadalta(l)))/(2*cdzaE0) ; 

coaff2  ■  (mv*«0''Cdz*v0*2*El+ind*dalta(l))/«2(i) ; 

poly2  ■  C(xgb(i)“2+zgb“2)  (-2*coaff2*zgb)  (coaff2“2-xgb(i)'2)] ; 

thataO  ■  asin(root8(poly2)); 

chaekl  ■  xgb(i)4cos(thata0(l))4zgb«sin(thata0(l)) ; 
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ch«ek2  «  xgb(i)<«>eoa(th«ta0(2))-*'Sgb*ila(tli«t«0(2)) ; 
if  •bs(eb*ek2-eo«ff2)  <  •b«(eb«ekl-eo«ff2) 
th«taer(i)  ■  tb«ta0(2); 

th«taer(i)  •  tbataOCl); 
and 

•Isail  daltad)  <  0 

wO  ■  -(zv*'sqrt(zv"2+4*edz*E04izd>t>d«lta(l)))/(-2*edz*E0) : 
coaff2  •  Cmw*w0<»^ed2««0*2*El-HBd*dalta(l))/«2(i) ; 

PQI72  ■  C(xgb(i)“2+*gb‘2)  (-2*eoaff2*zgb)  (co#ff2“2-xgb(i)‘2)] ; 
thataO  •  aBin(roots(pol72)) : 

ehaekl  ■  xgb(i)*eos(thataO(l))-fzgb«ain(thataO(l)} ; 
ehaek2  ■  xgb(i)*eoa(tbataO(2))-<>zgb*sia(th#taO(2)) ; 
if  ab8(ehack2-coafl2)  <  aba(ehaekl-coall2) 
thatacr(i)  •  tbata0(2); 

alza 

thataer(l)  •  thataO(l); 
and 
and 

and 

rasultz  >  Deg*  neritC:,!)  «0*onas(langth(xg) ,1)  thataer'] 

aava  aezit0.dat  xasnlts  -ascii 
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c 

c 

c 

c 


c 


c 

c 


PROGRAM  HOPP.O.POR 

EVALUATIOH  OF  HOPF  BIFURCATIOH  FORMDLAS 
ZERO  DIVE  PLAHE  ANGLE 


IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  L.n.MASS.MQDOT.MWDOT.NOl.TBETA. 


1 

2 

3 

4 

5 

1 

2 

3 

4 

6 
6 
7 


MQ .MV.MD .MDS .MDB .K1 .K2 . 

ALFA .BETA .GAMA .DELTA , 
B0.E1.E2.E3.E4. 

DW1.DU2.DV3.DW4. 

DQl.Dq2.Dq3.DQ4 

DOUBLE  PRECISION  M11.M12.M13.M21.M22.M23. 
M31.M32.M33. 

Nil .N12 .N13 .N21 .N22 .N23 . 
N31.N32.N33. 

L21 . L22 . L23 . L24 . L31 . L32 . L33 . L34 
L26 .L26 .L27 .L36 .L36 .L37 . 

L21A .L22A.L23A.L24A .L31A . 
L32A,L33A.L34A 


9 


DIMENSION  A(3,3)  .T(3,3)  ,TINV(3.3)  .FV1(3)  .m(3)  .TYT(3,3) 
DIMENSION  MR(3) ,WI(3) .TSATE(3.3) ,TLUD(3,3) ,IVLUD(3) ,SVLUD(3) 
DIMENSION  ASAVE(3.3) .Al(3.3) .A2(3.3) .ZL(26) .BR(26) 

DIMENSION  VEC0(25} .TEC1(26) .TEC2(26) .TEC3(25) .TEC4(26) 


OPEN  (20, FILE-* HOPF. RES '.STATUS- 'NEW*) 


WEIGHT-16S6.2363 


IT 

-  661.32 

L 

-  13.9792 

RHO 

-  1.94 

WRITE 

(*,*)  '  ENTER 

READ 

(♦,♦)  CD 

CD 

-  0.6*CD*RH0 

C 

-  32.2 

ZB 

-  0.0 

WRITE 

(♦,♦)  '  ENTER 

READ 

(*,•)  ZG 

ZG 

-ZG*L 

ZB 

-  0.0 

NASS 

-  WEIGHT/G 
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c 

c 

c 

c 


BOT  >  HEI6BT 
n>l  •  0.6*IIHO*L**2 
ZQOOT  •-6.3300E-04*0.6«RHO*L«*4 
ZWDOT  1 . 4629E-02*0 . 6*RH0*L**3 

ZQ  >  7.6460E-03*0.6*RI0«L**3 
ZU  >- 1 . 3910E-02«0 . B*RB0«L«*2 
MQOOT  «-8.8000E-04«0.6«RH0*L**5 
HWDOT  ■>6.6100E-04«O.S*REO*L*4>4 
HQ  *-3.7020E-03*0.6*RH0*L**4 
HH  -  1.0324E-02*0.6*RB0*L**3 

Z6B-ZG-ZB 

DEFIHE  THE  LEHGTH  X  AMD  BREADTH  B  TERNS  FOR  THE  IMTEGRATIOB 


IL(  D- 

0.0000 

IL(  2)- 

0.1000 

IL(  3)- 

0.2000 

XL(  4)- 

0.3000 

IL(  5)- 

0.4000 

XL(  8)- 

0.6000 

IL(  7)- 

0.6000 

XL(  8)- 

0.7000 

IL(  9)- 

1.0000 

ZL(IO)- 

2.0000 

IL(ll)- 

3.0000 

IL(12)- 

4.0000 

XL(13)« 

7.7143 

IL(14)- 

10.0000 

ILCIS)- 

15.1429 

IL(16)- 

18.0000 

XL(17)- 

17.0000 

IL(18)- 

18.0000 

ILC19)- 

19.0000 

IL(20)- 

20.0000 

IL(21)- 

20.1000 

IL(22)« 

20.2000 

IL(23)- 

20.3000 

IL(24)» 

20.4000 

IL(26)- 

20.4167 

DO  102  N  «  1,25 
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XL(II)  -  a/20.)*IL(M)  -  L/2. 


COirHUE 

BR(  D- 

0.000 

BR(  2)- 

0.486 

BR(  3)- 

0.668 

BR(  4)- 

0.778 

BR(  6)- 

0.871 

BR(  6)- 

0.946 

BR(  7)- 

1.010 

BR(  8)- 

1.060 

BR(  9)- 

1.180 

BR(IO)- 

1.410 

BR(ll)- 

1.670 

BR(12)- 

1.660 

BR(13)- 

1,670 

BR(14)- 

1.670 

BR(1S)« 

1.670 

BR(16)> 

1.630 

BR(17)» 

1.370 

BR(18)- 

0.919 

BR(19)« 

0.448 

BR(20)« 

0.196 

BR(21}« 

0.188 

BR(22)- 

0.168 

BR(23)- 

0.132 

BR(24)- 

0.063 

BR(2S)- 

0.000 

00  104  H 

-  1,26 

TECO(K)-BR(K) 

VEC1(K)-XL(K)*BR(K) 
TSC2(K)-ZL(K)*XL(K)*BR(K) 
VKC3(K)-XL(K)*XL(K)*XL(K)*BR(K) 
WC4(K)-XL(K)*XL(K)*XL(K)*XL(K)*BR(K) 
104  CORTIHUE 

CALL  IRAP(2S.VEC0.XL.E0) 

CALL  TRAP(25,TEC1,XL,E1) 

CALL  TRAP(25.VEC2.XL.E2} 

CALL  TRAP(26,VEC3,XL,E3) 

CALL  TRAP(2S,VEC4,XL,E4) 
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P^»**^SPP 


WRITE  (*.«)  '  ERTER  QAKKA' 

READ  (*.«)  EPRILOl 
Kaill— 0.01 
ZGIIAZ»0.01 

n6«80 

ZGKII-XQMIV*L 

Z6ifAX-Z6NAZ*L 

00  1  IT-l.IZG 
WRITE  (*,3001)  IT.IZG 
IG-ZGNIH*(ZGMiZ-ZGMII)*(IT-l)/(IZG-l) 

ZGB-ZG-ZB 

DT-(MASS-ZWOOT)*(IT-HqDOT)-(lfiSS*ZG*ZqOOT)*(MASS*ZG*NWDOT) 
CD6*C0/ ( 3 . 00*EPSILOM*DT) 

0W1-C06* ( ( IT-MqOOT) * ( -EO)*(llASS*ZG*ZqDOT) *E1 ) 

0W2-C06* ( ( n-MqOOT) * (3*E1) - (MASS*ZG*ZqDOT) *3*E2) 

0W3-C06* ( ( IT-NqOOT) * ( •3*B2) * (MASS*ZG*ZqDOT) *3*E3 ) 

DW4-C06* ( ( IT-NqOOT) * (E3 ) - (MASS*XG*ZqDOT) *E4 ) 

0qi-C06*((MASS-ZWD0T)*(El)*(HiSS*XG4NWD0T)*(*E0)) 

Dq2«C06*((NiSS>ZW00T)*(-3*E2)*OUSS*XG*mn>0T)*(3*El)) 

Dq3*C06*((MASS-ZWD0T)*(3*E3)*OfiSS*XG*HW00T)*(-3*E2)) 

Dq4«C08*((MASS-ZW00T)*(-E4)*(HAS8*XG4MWD0T)*(E3)) 

1HET10*ATAI(-XGB/ZGB) 

UO-(NASS-ZWDOT)  *  (n-NqOOT)  >OlWD0T*lfASS*XG)  *(ZqOOT*lfASS*XG) 
BBO- ( -ZWDOT*llASS-NASS*NW-Zq*MASS ) *XG 
*  *(-MASS*Nq*ZWDOT*Nq-ZqDOT*MW 

ft  -Zq*lfWDOT-HASS*NWDOT-IT*ZW*NqDOT*ZW) 

OCO— HASS*ZW*XG*Hq*ZW-Zq*HW-MASS*MW 
CC1-(>HASS*XG*ZW00T*XG*NASS*XB-ZWD0T*XB) *SI1 (TBETAO) 
ft  *( -NASS*ZB-ZWDOT*ZG*ZWDOT*ZB*MASS*ZG) *COS (TBETAO) 

0D1-(ZW*XG-ZW*XB)*SIB(THETA0)*(ZW*ZB-ZW*ZG) 
ft  *COS(THETAO) 

C 

WEI-BBO*CCO/ (AA0*DD1-BB0*CC1) 

TPO-DSqRT(  1566 . 2363/WEI) 
t^UO 
C 

C  DEIERNIHE  [A]  AND  CB]  COEFFICIEBTS 

C 

AllDT-(IT-MqOOT) *ZW*(MASS*X6*ZqD0T) *MW 

A12DT- ( IT-MqOOT) * (NASS*Zq) * (NASS*X6+ZqD0T) * (Hq-MASS*X6) 

A13DT—  (NASS*X6*ZqD0T)  *WEIGBT 
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i210T-(lflSS-ZVDOT)*MV«(NASS*X64mn>OT)  *ZU 
i22DT-(IUSS-ZVD0T)*(Mq'MASS*X6)’i>(MiSS*X64lIHl>OT)*(IUSS-i^Zq) 
A23DT—  (MASS-ZV00T)*ffEI6ET 
C 

ill-AllOT/DT 

A12-A120f/0V 

A13-A130V/DT 

A21«A2iOV/DT 

A22-A220V/DV 

A23-A230T/0V 

C 

Cl  lOT-  ( n-NqooT)  *mass«zg 
C120T— (HASS*X64ZqOOT) *MASS*Z6 
C210V—  (MASS-ZUDOT)  «NASS*ZG 
C220V-(MASS*XG4>MHD0T)  *MASS*ZG 
C 

Cll-CllDT/DV 

C12-C12DT/DV 

C21-C21DT/DT 

C2>C220T/0T 

C  EVALUATE  TRAXSFOBNATIOE  KATUX  OP  EIGEHTECTOES 

C 

El—  (XGB*SIH(THETAO)  -Z(ai*COS(TBZTAO)  ) 

K2— (1 . /6 . )*(ZGB*COS(THETAO) -XGB*SIH(THETAO) ) 

C 

A(l,l)-0.0 
A(l,2)-0.0 
A(l,3)>1.0 
A(2.1)-A13*K1 
A(2.2)-A11*U 
A(2.3)-A12*U 
A(3,1)-A23*K1 
A(3.2)-A214>U 
A(3,3)-A22*U 
DO  11  1-1,3 
DO  12  J-1.3 

ASATE(I.J)-A(I,J) 

12  CONTINUE 
11  CONTINUE 

CALL  RG(3,3,A,VR,VI,1,TTT,IT1,FT1.IEIUI) 

CAU  DSOMEGdEV.WB.WI, OMEGA, CHECK) 
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ONlGAO-OlflQi 
DO  6  1-1,3 

Td.D-  TTTd.IE?) 

Td,2)— TTTd.m+1) 

6  CQITim 

ZF  (XEf.EQ.l)  60  TD  13 
IF  (IST.EQ.2)  GO  TO  14 
STOP  3004 
14  DO  6  1-1,3 

Td,3)-TTTd,l) 

6  COITIIOB 
GO  TO  17 
13  DO  16  1-1,3 

Td,3)-md,3) 

16  COITIIUE 

17  COITUDE 

■OBNALIZATIOV  OF  THE  CUTICAL  HOSITKCTOR 

nOBM-l 

ZF  dlOBM.KE.O)  CiU  KORKALd) 

ZVYERT  IRAISFOBMATIOI  MATRIX 

DO  2  1-1,3 
DO  3  J«l,3 
TIIT(I,J)-0.0 
TSA?Ed,J)-Td,J) 

3  COMTIHUE 
2  COITIlUB 

CALL  DLin)(3,3,TSATE,3.TL00,I?LTJD) 

DO  4  Z-1.3 

ZF  dYLODCD.EQ.O)  STOP  3003 

4  COITIIUE 

CALL  0ILU(3,3.TLUD,ZTLUD,STLUD) 

DO  8  1-1,3 
DO  9  J-1,3 
TII?d,J)-TLUDd,J) 

9  COITIIUE 

8  COITIIUE 


CHECK  i:av(T)*A*T 


c 

nfOLT-l 

IF  (IHDLT.iq.l)  CALL  KDLTv:iMT,iSATB,T,A2) 

IF  (INULT.EQ.O)  STOP 

P-A2(3.3) 

PEIG-P 

C 

C  OEFIiniOH  OF  lij 

C 

Mll-TIHFd.l) 

II31-TIHT(2.1) 

iai-TIIY(3.1) 

H12-TIHV(1.2) 

II22«TIHT(2.2) 

II32«TIIF(3,2) 

M13-TII?(1.3) 

>23-TIST(2.3) 

H33-TIHV(3,3) 

C 

C  SEFUmOl  OF  Mlj 

C 

Mll-T(l.l) 

M21«T(2.i) 

M31-T(3.1) 

IU2«T(1,2) 

IQ2«T(2.2) 

M32«T(3.2) 

M13-T(1.3) 

If23-T(2,3) 

M33-T(3,3) 

C 

C  OEFINITIOH  OF  Lij 

C 

L26«C11*M31*N314C12«N21*N31 

L26>2*Cll*lf31*M32-K:i2*(If21*N324M22*H31) 

L27-C11*M32*H324C124'H22*M33 

L35-C22*N31*H31+C21*N21>»K3i 

L36>2*C22*H31*M32't-C21*(H21<t>H32-i-lf22*M31) 

L37>C22*N32*M32-^C21*H33*H22 

C 

C  DEFIHITIOH  OF  ALFA,  BETA.  6ANA 

C 
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D1  *1132*126  *  II33*L3S 
02  ->S2*L26  *  I3S*L36 
03  i4r32*L27  *  133*137 
C 

Oil— P 
012-0MBQA0 
021— 2*0MB6A0 
022— P 

023*2*0116010 
032*-0MB6A0 
033— P 
C 

BBTA*(02-021*01/011-023*D3/033)/(022-021*012/011-023*032/033) 
ALFA*(01-012*Bm)/011 
aAIIA*(03-032*BBTA)  /033 
C 

L21A*2*C11*ALFA*M31*M33*C12*ALFA*(M21*H334M23*M31) 

C 

L22A*2*C11*ALFA*M32*I133  *  2*C11*BBTA*M31*M33 
A  *  C12*ALFA*(I(22*I1334H32*II23) 

A  4  C12*BETA*(ll21*l»34N23*N3t) 

C 

L23A*2*Cll*0ANA*H31*113342*Cll*Bm*M32*M33 
A  4  C12*0ANA*(N21*N334||23*II31) 

A  4  C12*Bm*(ll22*N334N23*lf32) 

C 

L24A*2*Cll*0AMA*H32*N334C12*0ANA*(lf22*N334H23*N32) 

C 

L31A*2*C22*ALFA*N31*N334C21*ALFA*(M21*N334M23*H31) 

C 

L32A*2*C22*ALFA*M32*K3342*C22*BKTA*lf31*M33 
A  4  C21*ALFA*(K22*H334H32*N23) 

A  4  C21*BETA*(K21*M334M23*H31) 

C 

L33A*2*C22*6AIlA*lf31*N3342*C22*BBTA*M32*N33 
A  4  C21*0ANA*(H21*N334M23*N31) 

A  4  C21*BBTA*(H22*M334N23*II32) 

C 

L34A*2*C22*0AMA*N32*K334C21*6ANA*(M22*H334M23*N32) 

C 

L21*L21A4A13*K2*Mll**340«l*lf21**340H2*M31*M21**2 
A  4  1W3*II21*H31**240H4*M31**3 


n 


L22-L22A-»3*il3*K2*M12*Mll**2«S*0Hl«lf22*M21**2 
ft  *  Dil2*(2*l(21*l(22*K3l4H32*M21**2) 

ft  0H3*(2*il21*N31*H32+ll22«K31**2) 

ft  *  3*0W4*N32*M31*«2 

L23-L23A-i^3*A13*K2*Mll*M12«*2<»3*0Hl*M21*N22*4i2 
ft  IMf2*(M31*M22*«2«2*M21*H22*i(32) 

ft  0U3*(M21*M32**2«2*M22*lf31*M32) 

ft  *  3«0V4*H31*N32**2 

L24-L24A-»Ai3*K2*M12*«3-»0Hl*M22**3«DH2*ll32*if22**2 
ft  0H3*M22*M32**2+Dlf4*H32**3 

L31>L31A+A23*K2*Mll**3<»0Ql*lt21**3<»Dq2*M31*M21**2 
ft  *  0q3*M21*M31*«2-»l)Q4*N31**3 

L32«L32A+3«A23*K2*N12*N11**2+S*DQ1*M22*I121**2 
ft  *  Dq2*(2*H21*N22*lf31+lt32*H21**2) 

ft  ♦  I)q3*(24iM21*M3i*M32+lt22*K31**2) 

ft  *  >DQ4*N32*M31**2 

L33-L33A-«>3*A23*K2*Nll*M12**2-f3*DQl*ll21*M22**2 
ft  *  0q2*(M31*M22**2-»2*ll21*M22*ll32) 

ft  *  0Q3*(I121*M32**2-»2*lf22*H31*ll32) 

ft  •«>  3*0q4*N31*lfS2**2 

L34-L34A-»A23*K2*M12**34Dqi«K22**3^DQ2*II32*M22**2 
ft  *  0q3«N22*M32**24^Dq4*M32*«3 

C 

R11*M12*L21-M13*L31 

R12-H12*L224|13*L32 

ai3-H12*L23+N13*L33 

R14-M12*L244>H13*L34 

B21-N22*L21+M23*L31 

II22«N22*L22>N23*L32 

R23-N22*L23<^R23*L33 

II24-II22*L24«'I23*L34 

C 

C  EVALUATE  DALPHA  AHO  OONEGA 

C 

UIHC-0.001 
UR  -U-»UIHC 
UL  -U-UIHC 
U  -UR 
C 

A(l,l)-0.0 

A(l,2)-0.0 

A(1.3)-1.0 
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i(2.i)»il3*Ki 

i(2.2)-All*a 

A(2,3)«il2*U 

A(3.1)-A23*K1 

A(3.2)>A21*U 

A(3.3)-A22*0 

CAU  HC(3.3,A,HR,tfI.0.TTT,m.F¥l,inUl) 

CAU.  DSTABL(DE0S,H1.HI.FBEQ) 

ALPHR-OEOS 

OHEGR-FEEQ 

0-UL 

A(l.l)-0.0 

A(1.2)-0.0 

A(1.3)-1.0 

A(2.1)-A13*K1 

A(2.2)-All*a 

A(2.3)«A12*U 

A(3.1)-A23*K1 

A(3.2)«A21*U 

A(3,3)«A22«U 

CALL  B6(3,3,A,W»,WI.0,TTT.ITl,m,IEBR) 

CAU  DSTABLCOEOS.WR.HI.FBEq) 

ALPHL-DEOS 

OHEGL-FREq 

OALPSA-(ALPilR-ALPBL)  /  (UR-UL) 
OaMEGA-CONEOR-OMEGL) / (CR-DL) 

EFALUATIOS  OF  HOPF  BIFURCATIQR  CQEFFICIEITS 

COEFi-3 . 0*Rll+R13-»R22'i'3 . 0*R24 
C0EF2>3 .0*R21+R23-R12-3 .0*Ri4 
AMU2  — C0EF1/(8.0«0ALPHA) 

BETA2-0.2S*CQEFi 

TATI2  — (C0EF2-0QlfEGA«C0RFl/DALPHA}/(8 .0*0HE6A0) 
PER  >2 . 0*3 . 1416927/OHEQAO 
PER  -PER*U/L 

WRITE  (20.2001)  Z6/L.(J.COEF1,DALPHA.ONE6AO.PEI6 


1  coirnuE 
STOP 

1001  FOUUT  ('  miR  lOMBEI  OF  DATA  LUES*) 

1002  FOSMAT  (*  EHTEE  00.  26.  AMO  DSAT') 

1003  FOIMAT  (>  EHTEE  BOH  PLAIE  TO  STEIl  PLAHE  EATIO'} 

1004  FOEMAT  (*  EHTEE  ZQ') 

2001  FOEMAT  (6E14.6) 

4001  FOEMAT  (3F16.6) 

3001  FOEMAT  (2IS) 

EHD 

SOBEODTIHE  DS0ME6(IJX.VE.WI. OMEGA. CHECK) 

IMPLICIT  DOUBLE  PEECISIOH  (A-H.O-Z) 

OIMEHSIOH  UE(3).HI(3) 

CHECK— 1.0042S 
DO  1  >1.3 

IF  (HECD.LT.CHECK)  GO  TO  1 
CHKCK-HE(I) 

U-I 

1  COHTIHDE 
OMKGA>DABS(VI(IJ)) 

IF  (HI(IJ).6T.0.00)  UK-IJ 
IF  (HI(IJ).LT.O.OO)  UK-IJ-1 
BETUBI 
EHD 


SOBEODTIHE  DSTABL(OEOS.VE.HI. OMEGA) 
IMPLICIT  DOUBLE  PEECISIOH  (A-H.O-Z) 
OIMEHSIOH  UE(3).HI(3) 

OEOS— 1.00<»20 
DO  1  >1.3 

IF  (HE(I).LT.DEOS)  GO  TO  1 
OEOS-WE(I) 

U-I 

1  COHTIHUE 
OMEGA-WI(IJ) 

OMEGA-DABS(OMEGA) 


EHD 

CassaaaBaaaaassaaaaaaaaaBaaaatBaaasaaBaBBmiaBaBsaBBaasBBBBBBBaBsa 

SOBEODTIHE  HOEMAL(T) 

IMPLICIT  DOUBLE  PEECISIOH  (A-H.O-Z) 


DOiBISIOH  T(3.3).TI0a(3,3) 

CrAC»T( 1 . 1 ) **2+T( 1 , 2) **2 

IF  (DABSCCFAO.LE.Cl.D-lO))  STOP  4001 

THOR(1.1)>1.DO 

T10R(2,1)«(T(1,1)*T(2.1)+T(2.2)*T(1.2))/CFAC 

TI0R(3.1)-(T(1.1)*T(3.1)+T(3.2)*T(1.2))/CFAC 

‘niQR(l,2)-0.00 

TR0R(2.2)-(T(2,2)*T(1,1)-T(2.1)*T(1.2))/CFAC 
TI0R(3.2)-(T(3,2)*T(1.1)-T(3,1)*T(1,2))/CFAC 
DO  1  >1.3 
00  2  J-1.2 

T(I.J)»THOR(I,J) 

2  COITIVOB 

1  COITIRUE 
BSTORH 
EHD 

SDBROUTIHB  HDLT(TIHV.A.T.A2) 

IMPLICIT  DOUBLE  PRECISIOI  (A-H.O>Z) 

DIMBISIOH  T1ET(3.3),A(3,3),T(3,3),A1(3,3),A2(3,3) 
DO  1  >1,3 
DO  2  >1,3 
A1(I,J)«0.D0 
A2(I,J)-0.D0 

2  COITIIIUE 
1  CORTHUE 

DO  3  >1,3 
DO  4  >1,3 
DO  6  >1,3 

A1(I.J)-A(I.K)*T(K,J)+A1(I,J) 

5  COITIIIUE 
4  COITIHUE 

3  COITIHUE 
DO  6  >1.3 

DO  7  >1,3 
DO  8  K-1.3 

A2(I,J)-TII7(I,K)*A1(K.J)+A2(I.J) 

8  COITIHUE 

7  COITIHUE 

6  COHTINUE 

DO  11  >1.3 

C  WRITE  (•,101)  (A(I,J),>1,3) 
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11  COMTIHUB 

DO  12  1-1.3 

C  MMTE  (*.101)  (T(I,J).J-1,3) 

12  COMTIMUE 

00  10  1-1.3 

C  WRITE  (*.101)  (i2(I.J).J-1.3) 

10  coirnuE 

C  WRITE  (*.101)  12(1.1) 

RSTURI 

101  FORMIT  (3B16.6) 

ERD 

SOBROOTIHE  TRAP(H.A.B.OnT) 

C 

C  NOMERICIL  INTE6R1T10H  ROUTIME  USIIG  THE  TRAPEZOIDAL  ROLE 
C 

IMPLICIT  DOUBLE  PRECISIOI  (A-H.O-Z) 

OIMEHSIOH  A(1).B(1) 

Hl-M«l 
OOT-0.0 
DO  1  I-l.Nl 

OOT1*0,6*(A(I)+A(I+1))*(B(I+1)-B(I)) 

OUT  -OUT+OUTl 

1  comthue 

RBTURH 

EHD 


PROQlUlf  BOPF.DELTA.FOR 

ETALOATXOH  OF  HOPF  BIFURCATIOII  FORMULAS 

STBRR  PLANE  ANGLE  IS  NON  ZERO 


IMPLICIT  DOUBLE  PRECISION  (A-B.O-Z) 

DOUBLE  PRECISION  L. IT. MASS .MQDOT.MWDOT.NDl, THETA, 


1 

2 

3 

4 
6 

DOUBLE  PRECISION 

1 

2 

3 

4 
6 
« 
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MQ .MV ,MD .ZD .MDS ,MDB .K1 .K2 , 

ALFA .BETA .GAMA .DELTA . 
E0.E1.E2.E3,E4. 

DV1.DW2.DV3.DV4. 

Dqi.DQ2.DQ3.DQ4 

M11.M12.M13.M21.M22.M23, 

M31.N32.M33. 

Nil . N12 .N13 .N21 .M22 .N23 , 
N31.N32.N33. 

L21 . L22 .L23 , L24 . LSI . L32 . L33 . L34 . 
L26 . L26 , L27 . L3S . L36 . L37 . 

L21A . L22A .L23A . L24A . L31A , 
L32A.L33A.L34A 


DIMENSION  A(3,3),T(3.3),TINf(3.3),F¥l(3),m(3),TTT(3,3) 
DIMENSION  WR(3),VI(3).TSATE(3,3).TLUD(3.3).ITLUD(3).STLUD(3) 
DIMENSION  ASATE(3 .3) ,A1(3.3) .A2(3.3) .IL(25) .BR(2S) 

DIMENSION  VEC0(25).TEC1(26).7EC2(26),VEC3(26),7EC4(26) 


OPEN  (10. FILE- 'HOPF. DAT '.STATUS'* OLD >} 
OPEN  (20 , FILE- ' HOPF .RES', STATUS- ' NEV ' ) 


HEIOBT-1666.2363 

n  «  561.32 

L  -  13.9792 

RHO  -  1.94 

WRITE  (*.«)  '  ENTER  CD' 

READ  (*,*)  CD 

CD  -  0.6*CD*RH0 

0  -  32.2 

ZB  -  0.0 

WRITE  (*,•)  '  ENTER  ZG/L' 

READ  (*,*)  ZO 

ZG  -ZG*L 

ZB  -  0.0 
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HASS  «  WEI6BT/G 
BOT  -  WEIGHT 
MDl  -  0.5*RH04>L**2 
ZqOOT  «-6.3300E-04*O.S*RHO*L**4 
ZWDOT  — 1 . 4629E-02*0 . 6*RH0*L**3 
Zq  -  7.6460E-03*0.6*RB0*L**3 
ZW  — 1.3910E-02*0.6*RH0*L**2 
MqOOT  >-e.8000E>04*0.6*RHO*L**6 
HWDGT  — 6.6100E-04*O.S«ltHO*L**4 
Mq  — 3.7O2OE-O3*O.S*RH0«L**4 
HW  «  1.0324E-02*O.S4>RHO*L**3 
C 

ZGB-ZG-ZB 


C 

C  DEFINE  TIE  LENGTH  X  AND  BREADTH  B  TERMS  FOR  THE  INTEGRATION 
C 


XL(  D- 

0.0000 

IL(  2)- 

0.1000 

IL(  3)- 

0.2000 

IL(  4)- 

0.3000 

ZL(  6)« 

0.4000 

ZL(  6)> 

0.6000 

XL(  7)- 

0.6000 

IL(  8)- 

0.7000 

IL(  9)- 

1.0000 

IL(IO)- 

2.0000 

IL(11)« 

3.0000 

ZL(12)« 

4.0000 

ZL(13}- 

7.7143 

ZL(14)- 

10.0000 

XL(16)- 

16.1429 

ZL(16)- 

16.0000 

XL(17)» 

17.0000 

IL(18)- 

18.0000 

IL(19}- 

19.0000 

IL(20)- 

20.0000 

IL(21)> 

20.1000 

ZL(22)> 

20.2000 

ZL(23)- 

20.3000 

ZL(24)« 

20.4000 

IL(26)« 

20.4187 
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DO  102 

N  -  1.26 

XL(H) 

-  a/20.)*XL(H)  -  L/2 

CORTIHUE 

BR(  D- 

0.000 

BR(  2)- 

0.485 

BR(  3)- 

0.668 

BR(  4). 

0.778 

BR(  6)- 

0.871 

BR(  6)- 

0.946 

BR(  7)- 

1.010 

BR(  8)- 

1.080 

BR(  9)- 

1.180 

BR(10)> 

1.410 

BR(ll)- 

1.570 

BR(12)> 

1.860 

BR(13)- 

1.670 

BR(14)- 

1.670 

BR(16)- 

1.870 

BR(16)- 

1.630 

BR(17)- 

1.370 

BR(18)- 

0.919 

BR(19)« 

0.448 

BR(20)- 

0.196 

BR(21)» 

0.188 

BR(22)- 

0.168 

BR(23)- 

0.132 

BR(24)> 

0.053 

BR(2S)- 

0.000 

DQ  104  K 

-  1.25 

VECO(K}-BR(K) 

TBC1(K)-XL(K)«BR(K) 
VEC2(K)-XL(K)*XL(K)*BR(K) 
VEC3(K)-XL(K)*XL(K)*IL(K)*BR(K) 
VRC4(K)-XL(K)*XL(K)*XL(K)*XL(K)*BR(K) 
104  CORTmUE 

CALL  TRAP(25.VEC0,XL.E0) 

CALL  TRAP(26,VEC1,XL.E1) 

CAU  TRAP(26,VEC2.XL.E2) 

CALL  TRAP(26,VEC3,XL,E3) 
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CALL  T11AP(26.VEC4.XL.E4) 

WRITE  (*.*)  '  ENTER  GAHMA* 

READ  (*.*)  EPSILON 
IXG-80 

CaaaaaaBaaMaaBaaaaaaaaaaMaaaaauaBBuamaaBBBraBBaaBMsaaaaBssBs 

00  1  n-i.ixG 

READ  (10.*)  XG.U.WO.THETAO 
WRITE  (*,3001)  IT.IXG 
X6«XG*L 
XGB-XG-XB 

0T>(MASS-ZWD0T) * (n-HQOOT) -(MASS*XG*ZQDOT) * (MASS*XG*MWDOT) 
C0NST2-CD/ (3 . 0*0T*EPSIL0N) 

TDW1-C0NST2*( (IT-MQDOT)*(-EO)+(MASS*XO+ZQDOT)*(E1) ) 
TDW2-C0HST2* ( ( IT-MQDOT) * (3 . 0*E1 ) +(MASS*XO+ZQDOT) * ( -3 . 0*E2 ) ) 
TDW3-C0NST2* ( ( IT-HQDOT) * (-3 . 0*E2) * (HASS*XG*ZQDOT) * ( 3 . 0*E3 ) ) 
TDW4-C0NST2* ( ( IT-MQOOT) * (E3) *(MASS*X6*ZQD0T) * ( -E4) ) 
TDqi-C0HST2*((MASS-ZWD0T)*(El)*(NASS*XG*MWD0T)*(-E0)) 
TDq2-C0HST2* ( (MASS-ZWDOT)* ( -3 . 0*E2) * (MASS*X0*NWD0T) * (3 . 0*E1 ) ) 
TDq3-C0NST2* ( (MASS-ZWD0T)*(3 . 0*E3) 4(MASS*XO*MWDOT) * (-3 . 0*E2) ) 
TDq4*C0NST2*( (HASS-ZWD0T)*(-E4)*(HASS*XG4MWD0T)*(E3) ) 

C 

CONST«TANH(WO/EPSILON)*CD/DV 

0W1«C0HST*(  (IT-MqOOT)*('EO)  *(MASS*X6*ZqD0T)*El) 

0W2-C0NST*(  (IT-Mq00T)*(2*El)  >OlASS*X6*Zq00T)*2*E2) 

DW3«C0NST*(  (IT-MqD0T)*('E2)  *(NASS*XG*ZqD0T)*E3) 

Dqi-CONST*(  (HASS-ZW00T)*(E1)  HNASS*XG*MWDOT)*(-EO)) 

Dq2-C0NST*(  (HASS-ZWD0T)*(-2*E2)+(HASS*X6*NWD0T)*2*Ei) 
0q3-C0NST*(  (HASS-ZW00T)*(E3)  *(NASS*XG*MWD0T)*(-E2)) 

C 

C  DETERMINE  CA]  AND  [B]  COEFFICIERTS 

C 

A1 lOT- ( IT-MqOOT) * ( ZW-2*CD*E0*W0) 

A  *(NASS*XG*ZqOOT)*(lfW*2*CD*El*WO) 

A12DT-(IT-MqOOT)*(MASS4Zq*2*CD*El*WO)*(MASS*Xa*ZqOOT) 
k  *(Mq-NASS*XG-HASS*ZG*W0-2*CD*E2*N0) 

A130V— (HASS*XG*ZqDOT) *WEIGHT 
A21DT-(MASS-ZWDOT)*(HW*2*CO*E1*WO) 

A  *(NASS*XG*MWD0T)*(ZW-2*CD*E0*W0) 

A220T-(MASS-ZWDOT)*(Nq-HASS*XG-HASS*ZG*WO->2*CD*E2*WO) 

A  +(HASS*XG*NWDOT)*(MASS*Zq42*CD*El*WO) 

A23DT— (MASS-ZWOOT) *WEIGET 
C 
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jpj  Ji.WUi 


ill-AllDT/OV 
A12-A120T/DV 
A13-A130T/0V 
A21-A21DV/DT 
A22-A220V/DV 
A23-A23DT/DT 
C 

CllDV- ( IT-HQDOT) «MASS«ZG 
C12DV— (NASS*Za«ZQOOT) *MASS*ZG 
C21DV— (NASS-ZVDOT) *MASS*ZG 
C22DT-  (NASS*ZG4MWOOT)  *NASS4>ZG 
C 

Cll-CllOV/DV 
C12-C120T/0T 
C21-C21DT/DV 
C22-C220V/0V 

C  E7ALUATE  TRAISFORMATIOH  MATRIX  OF  EIGBITECTORS 

C 

ICl«-(IGB*SII(TlIETAO)-ZGB*COS(THETAO)  ) 

K2—(l . /6 .  )*(ZGB*COS(TBETAO)-IGB*SIl(THETAO)) 

C 

A(l,l)-0.0 
A(1.2)-0.0 
A(1.3)-1.0 
A(2.1)-A13*K1 
A(2.2)-A11*U 
A(2.3)-A12*U 
A(3.1)-A23*K1 
A(3.2)-A21*V 
A(3.3)-A22*U 
DO  11  1-1,3 
DO  12  J-1,3 

ASAVE(I.J)-A(I.J) 

12  OOMTIlUE 
11  COITIHUE 

CALL  RG(3.3.A.HR,UI,1.TTT.IT1,FT1.IXRR} 

CALL  DSOMEGCIEV.VR.WI, OMEGA, CHECK) 

OMEGAO-OMEGA 
DO  6  >1,3 
T(I,1)-  TTT(I,IE?) 

T(I,2)— TTT(I,IEV+1) 
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coimuE 

IF  (lET.EQ.l)  GO  TO  13 
IF  (lET.EQ.S)  GO  TO  14 
STOP  3004 
DO  6  1-1.3 

T(I.3)-TTT(I.l) 

COITIIUE 
GO  TO  17 
00  16  1-1.3 
T(I.3)-TTT(I.3) 

COITIHUE 

COETIMUE 

HOlWALIZATIQE  OF  THE  CRITICAL  EKXITECTOR 
nORM-l 

IF  (IIORM.NE.O)  CALL  RORIfALCT) 

hfert  traisfqrmatioh  katrix 

DO  2  1-1.3 
DO  3  J-1.3 
TIIY(I.J)-0.0 
T5AVE(I.J)-T(I.J) 

COITIiniE 

coithoe 

CALL  DLUD(3.3.TSATE.3.TL0D.IVLUD) 

DO  4  1-1,3 

IF  (ITLUO(I).Eq.O)  STOP  3003 
COITIIUE 

CAU  DILU(3.3,TLUD.ITLUD.SFLUD) 

DO  8  1-1,3 
DO  9  J-1,3 
TIIT(I.J)-TLUD(I.J) 

COITIHUE 

COITIIUE 

CHECK  Ixit(T)*A*T 
DfULT-l 

IF  (IMULT.Eq.l)  CALL  NDLT(TIIT.ASATE.T.A2) 
IF  (IMULT.EQ.O)  STOP 


IP-A2(3.3) 

PBI6-P 

C 

C  KFIMITIOM  OP  llj 

C 

■ii«Tiinr(i,i) 

«3i-Tiinr(3,i) 

■12-TIMV(1,2) 

M22-TII¥(2,2) 

■3>Tm(3,2) 

>13-TIIV(1.3) 

I2S«TIIT(2.3) 

II33-TIIT(3.3) 

C 

C  OEPIMniOI  OF  Mlj 

C 

M21«T(2,i) 

K31-T(3,l) 

M12-T(l.2) 

M22-T(2,2) 

I132«T(3.2) 

M13-T(l,3) 

M23-T(2.3) 

IB3*T(3,3) 

C 

C  OEFlHmOI  OF  Uj 

C 

I.2S-Cll*N31*M3i4C12*N21*ll31 

ia6-2*Cll*M31*M32+C12*(II21*ll32+M22*M31) 

I37»C11*!I32*M32+C12*I!22*M33 

L36«C22*ll3i*lf31^C21*II21*ll31 

I36«2*C22*H31*II32>G21*  (ll2i4'll32«N22*lf31 ) 

LS7>C22*N32*M32<i’C21*H334i|Q2 

C 

C  DEFIVITIOI  OF  ALFA,  BETA,  GAHA 

C 

01  •N32*L26  *  lf33*L36 
02  «II32*L26  *  I33*L36 
03  -M32*L27  *  133*1.37 


C 
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Dll— P 
D12-0NE6A0 
D21— 2*0ME6A0 
D22— P 

023-2*0IIE6A0 
D33— 0HB6A0 
D33— P 

C 

BITA-(02-D21*Dl/011-D23*D3/D33)/(D22-D21*012/011-D23*D32/033) 

ALFA-(D1-D12*BBTA)/011 

OAHA- (03-D32«Bm) /033 

C 

L21A-2*  (Cl  140W3  )  *  ALFA*M31*M33-»  (C12<fOH2  }  *ALFA*  (M21  *N33-»N23«N3 1 } 
A  >2*0W1*M21*II23«ALFA 

C 

L22A-2*(C1140U3)*ALFA*II32*N33  *  2*(Cll<i^0U3)*Bm*M31*N33 
A  *  (C124DU2)*ALFA*(M22«lf33-i>H32*ll23) 

A  *  (C12-»DW2)*BBTA*(M21*M334M23*M31) 

A  •*-  2*0Vl*(H21*N23*Bm+N22*M23*ALFA) 

C 

L23A-2*(Cll<»^0V3)«aAMA*ll31*lt33-»2*(Cll-»DV3)*BBTA*N32*M33 
A  *  (C124DV2)<»0AlfA*(lf21*i(33‘«^H23*ll31} 

A  *  (C124DH2)*BBTA*(ll22*if33-Mf23*ll32) 

A  *  2*0tfl*(N21*lf23«aAKA-»M22*H23*Bm) 

C 

L24A-2«(Cll+0H3)*aAMA«M32<tiN334(C12-i^DH2)*QAHA*(M22*M33-^lf23*M32) 
A  ^2«DVl*N22*N23«aANA 

C 

L31A-2*(C22^0q3)*ALFA*N31*M33-»^(C21’^DQ2)*ALFA*(M21*H33'»H23*M31) 
A  •»2*Dqi*N21*l(23«ALFA 

C 

L32A-2*(C2240q3}*ALFA*N32*M33  *  XC22-^0q3)*BBTA*N31*M33 
A  *  (C21<»0q2)«ALFA*(N22*N33+M32*ll23) 

A  *  (C2l40q2)*BBTA*(ll21*II334|f23*lf31) 

A  4  2*0qi«(lf21*H23*Bm4H224|f23*ALPA) 

C 

L33A-24 (C2240q3) *6AMA*N31*M33424 (C224Dq3) *BBTA4M32*H33 
A  4  (C2l4Dq2)*GANA*(N2l4M334H234H31) 

A  4  (C2l4Dq2)*BETA4(H22*N334M23*N32) 

A  4  2*Dqi*(N21*H23*GA}(A4lf22*N234Bm) 

C 

L34A-2* (C2240q3) *GAIfA*N32*N334 (C2l4Dq2 ) *GAMA* (M22*M334H23*M32) 
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pipippi 


k  ♦2*0qi*M22«N33*aAMA 

C 

L21-L21A<»il3*K3*Mll**3 
L22«L22A-»3*A13*I2*K12«M11*«2 
L23>L23A*3*A13*K2*N1 l*Mi2*«2 
L24>L24A4A13*K2«N12**3 
L31-L31A-»A23*K2*N11*«3 
U2-L32A-i^3*A23*K2*M12«Nll*«2 
L33-L33A-»3*A23*K2*Nll*Mi2««2 
L34-L34A-»A23*K2*Ni2*«3 
C 

L21«L2i-»TDWl*N2i«*3-<'TDH2*if3l4>ll21**2-fTDV3*M2i*H31**2’<-TDV4*M31**3 
L32>L22>TDH1*3 . 0«N22*II21**2-*>TDH2*(2 . 0*ll2i*ll31*M22<i-N32*lt21*«2) 

A  ♦TDV3*(2.0*ll31*ll32*lf2l4M22*N31**2)<MDH4*3.0*M32*M31**2 

L23>L23«TDH1*3.0*M21*H22**2’^TDV2*(M31*M22**2+2.0*M21*N22*M32) 

A  ♦TDira*  (M21*lt32**2<»2 . 0*II31*M32*N22)  ♦TDV4*3 . 0*N31*M32««2 

L24-L244TDV1*K22«*34TDV2«I132*M22«*2-»TDU3*K22*M32«*24'TDU4*H32**3 
L31-L31-»TDQl*ll21**3«TDQ2*lf31*ll21**24TDQ3*N21*N31**2-«-T0Q4*N31**3 
L32«L32>TDqi*3 . 0*M22«lf21**2-^T0q2*(2 . 0*ll21*N31*lf22-»M32*IQl**2) 

A  ♦TDq3«(2. 0«ff31*lf32*ll21+K22*il3i**2)-^TDq4*3 . 0*M32*N31**2 

L33-L33<»TDqi*3 . 0*M21*lf22««2«TDq2*(lf3i*M22**2<»2 .0*lt21*M22*N32) 

A  ♦TDq3*  (N21*H32*«2>2 . 0<i>N31*M32*N22)  ••>TI>q4«3 . 0*M31*M32«*2 

L34-L34<»TDqi*N22**3<fTDq2*ll32*lf234>*2i'TDq3*H22*ll32'«>4>2-^TDq4*N32**3 

C 

B11«V124>L21+I13*L31 

R12-M12*L224M13*L32 

R13-M12*L23-»M13*L33 

1114-1I12*L244H13*L34 

R21-H22*L21+H23*L31 

B22«I22*L22+I23*L32 

R23-H22*L23-»I23*L33 

B34->H22«L24«I23«L34 

C 

C  EfALUATE  DALPHA  AID  DONEQA 

C 

unc-o.ooi 
OR  ^-^unc 
DL  ^-UIIC 
U  ^ 

C 

A(l.l)-0.0 

A(1.2)-0.0 
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a  n 


A(1.3)>1.0 

i(2.1)-A13*Kl 

A(2.2)«A11*U 

A(2.3)-A12*U 

A(3.i)-A23*K1 

A(3,2)-A21*!; 

A(3.3)>A22*U 

C 

CALL  Bfl(3,3,A,Wl,WI,0,TTT,IYl,m,ma) 
CAU  OSTABLCOEOS.HE.UI.FBSq) 

A::.P1R"0B0S 

QMBOR-PKiq 


A(1.1)«0.0 

A(l,2)«0.0 

A(l,3)«1.0 

A(2,1)«A13*X1 

A(2.2)«Ali*U 

A(2.3)«A12*U 

A(3,1)«A23«K1 

A(3,2)«A21*V 

A(3.3)-A224>U 

C 

CAU  M(3,3,A,Wl,WI.0,TTT,m.m,IERa) 

CAU  DSTABL(DEOS.WR.HI.FBSq) 

ALPEL-DEOS 

OHEGL-FREQ 

C 

OALPHA- ( ALPHR-ALPHL) / (UR-DL) 
001IE6A«(0MB6R-0I1E6L)  /  (UR-QL) 

EfALUATIOI  OF  BOPF  BIFORCATIOI  COBFFICIEITS 
C 

C0EFi«3 . 0*Rll<»R13-»R22-i>3 . 0*R24 
C0EF2-3 . 0«R21>R23-R12-3 . 0*R14 
AMD2  »C0EF1/(8.0*DALPHA) 

BETA2-O.26*C0EFl 

C  TAD2  — (C0EF2-D0HE6A*C0EF1/DALPHA)/(8.0*01IE6A0) 

C  PER  -2.0*3. 1416927/0ME6A0 

C  PER  -PER*U/L 
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iOlTE  (20.2001)  X6/L.n,C0EF1.0ALPHA.01fEGAO.PEZ6 
1  COVTnniE 
STOP 

1001  FORMAT  ('  ERTER  BOMBER  OP  DATA  UMES') 

1002  FORMAT  (’  EBTER  00.  Z6.  AMO  DSATO 

1003  FORMAT  (*  ERTER  BOV  PLARE  TO  STEER  PLARE  RATIO') 

1004  FORMAT  ('  ERTER  Z6') 

2001  FORMAT  (6E14.6) 

4001  FORMAT  (3F18.8) 

3001  FORMAT  (218) 

ERD 

SOBROOTIRE  0S0MEG(IJK.HR.WI.0ME6A.CBECK) 

IMPLICIT  DOOBLE  PRECISIQR  (A>H.0-Z) 

OIMERSIOR  WR(3).WI(3) 

CHECK— 1.004-28 
00  1  >1.3 

IP  (VR(I).LT. CHECK)  GO  TO  1 
CHECK-HR(I) 

II-I 

1  CORTIHOB 

aMKGA«OABS(HI(IJ)) 

IF  (HI(IJ).6T.0.D0)  IJK-IJ 
IF  (WI(IJ).LT.O.OO)  UK-IJ-1 
RETORR 
ERD 

Oaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

SOBROOTIRE  OSTABL(DEOS.VR.VI. OMEGA) 

IMPLICIT  OOOBLE  PRECISIQR  (A-H.O-Z) 

OIMERSIOR  VR(3).HI(3) 

ISOS— 1.0D-»20 
00  1  >1.3 

IF  (HR(I).LT.OEOS)  GO  TO  1 
OBOS-HR(I) 

IJ«I 

1  CORTIROE 
0ME6A>HI(IJ) 

OHEGA-OABS(ONEGA) 

RETORR 

ERD 

(^aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

SOBROOTIRB  BORMAL(T) 
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O' 


IMPLICIT  DOUBLE  PRECISIQM  (i-H.O-Z) 

OIMBMSIOH  T(3.3).TMQR(3.3) 

CP1C»T( 1,1) **2+T( 1 , 2) **2 

IF  (DABS(CFiC).LE. (1.0-10))  STOP  4001 

TVOR(1.1)>1.00 

TI0R(2,1)-(T(1,1)*T(2.1)+T(2.2)*T(1.2))/CFAC 

TI0R(3.1)-(T(1.1)*T(3.1)+T(3.2)*T(1.2))/CFAC 

TVOR(1.2)-O.DO 

T10R(2.2)-(T(2.2)*T(1.1)-T(2.1)*T(1,2))/CFAC 
T*0R(3.2)-(T(3.2)*T(1.1)-T(3.1)*T(1.2))/CFAC 
DO  1  1-1,3 
00  2  J-1,2 

T(I,J)-TMOR(I,J) 

2  COHTIIIUE 
1  COMTIHUE 
BETURI 
EHD 


SUBROOTIIE  MDLT(TINT,A,T,A2) 

IMPLICIT  DOUBLE  PRECISIOM  (A-B,0-Z) 

DIMEBSIOM  nBV(3,3) ,A(3,3) ,T(3,3) ,A1(3,3) ,A2(3,3) 
DO  1  1-1,3 
DO  2  >1.3 
Al(I.J)-0.00 
A2(I.J)-0.00 

2  COBTIRUE 
1  COimiUE 

DO  3  1-1,3 
DO  4  J-1,3 
DO  5  K-1,3 

Ai(I,J)-A(I.K)*T(K,J)+Al(I,J) 

5  COITIIUE 
4  COITIHUE 

3  COMTIHUE 
DO  6  1-1,3 

DO  7  J-1,3 
DO  8  K-1.3 

A2(I,J)-TIHT(I,K)*A1(K,J)+A2(I,J) 

8  COMTIHUE 

7  COMTIHUE 

6  COMTIHUE 

DO  11  1-1,3 
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C  HRTTE  (*.101)  (A(I.J).J-1.3) 

11  corrnoE 

00  12  I«l,3 

C  WETTE  (*.101)  (T(I,J),>1,3) 

12  CQITIHtlB 

00  10  >1.3 

C  UUTE  (*.101)  (A2(I.J).J-1.3) 

10  corrnniE 

C  miTE  (*.101)  A2(l,l) 

RETOSM 

101  FOBMAT  (3E16.6) 

EIO 

SDBRODTIIE  TRAPd.A.B.OOT) 

C 

C  HDHERICAL  IHTEGRATIOI  BOOTIHE  USII6  TEE  TSAPEZOIOAL  BOLE 
C 

IMPLICIT  OOOBLE  PBECISIOI  (A-B.0-1) 

OIMERSIOE  A(1),B(1) 

Ml-B-1 
OOT"0.0 
00  1  >1.H1 

00Tl-0.8*(A(I)*A(I+l))*(BOl)-B(I)) 

OUT  ■OUT+OUTl 
1  COITIHUE 
BETOBM 
ERD 
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PR06RAH  SIM. FOR 


C 

c 

C  SINUUTIOH  PROGRAM  OSIHG  A  FOURTH  ORDER  AOAMS-BASHFORTH 

C  IMTEGRATION  SCHEME 

C 

REAL  L.MW.MASS.IT.MQOOT.MVDOT.MQ. 

*  U.TIME,HDl.ND2.irD3.in)4. 

*  F1(4).F2(4).F3(4) 

C 

DIMEMSIOH  XL(26) .BR(2S) .VEC1(26) .TEC2(26) 

C 

OPEM  (10, FILE- 'SIM. DAT ’.STATUS-* OLD') 

OPEH  (20. FILE- 'SIM. RES'. STATUS- 'HEV') 

C 

C  EHTER  SPEEDS  AHD  TIME  DATA 
C 

READdO.*)  U.TSIM.DELT.IPRMT. 

«  ZG.ZO.THETA.U.q 

C 

C  GEOMETRIC  PROPERTIES  AMD  HTDRODTHAMIC  COEFFICIEHTS 
C 

PI  -4.0*ATAM(1.0) 

THETA  -THETA*PI/180 

RHO  -1.94 

CDZ  -0 . 50 

L  -13.9792 

HDl  «O.6*RB0*L**2 

HD2  -0.6*RB0*L**3 

IID3  -0.6*RH0*L**4 

in)4  -0.6*RHO*L**5 

HEIGHT-1666 . 2363/ (ND1*U**2) 

MASS  -1566. 2363/ (32. 2*MD2) 

It  -661.32/H04 

ZQDOT— 6.3300E-04 
ZHDOT— 1.4629E-02 
Zq  -  7.5460E-03 
ZH  — 1.3910E-02 
MqOOT— 8.8000E-04 
MHDOT— 5.6100E-04 
Mq  — 3.7020E-03 
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Ml 


-  1.0S24B-02 


C 

ZB  *0.0 

ZB  -0.0 

Z6B  -ZG-ZB 
ZQB  -ZG-ZB 

C 

C  DBTBBNIIE  Ci]  AID  [B]  COBFriCIZITS 

C 

D?  -(MISS-ZVDOT)  *  (IT-HQDOT)  -(l(iSS*ZG4ZQD0T)*  (lliSS*ZG^MWDOT) 

AllDV-(IT-MQOOT)*ZV4(|fA88«ZG4^ZQOOT)*IW 
il20f-(IT-MqD0T)*(MASS-^ZQ)«(IULSS*ZG4ZQD0T)*(MQ-MASS*ZG) 
AISDT—  <IIASS*ZG^ZqDOT)  *HIIGBT 
A210T-(lfiS8-ZirD0T)*IW«(lli88*ZG4|fin)0T)*ZW 
A220T-(NASS-Zin)0T)*(Mq-l(iSS*ZG)-«'(ltASS*ZG-»Min)OT)«(MASS4Zq) 
A23DV— (MAS8-ZVDGT)*liBIGIT 
C 

All-illOT/DT 

A12-A120Y/0V 

il3-A13DT/DT 

A31-A2iOT/DV 

A22-A22DT/0Y 

123-A23DY/DT 


C 

C  DKFUE  the  LEIGTH  Z  aid  BIEADTH  B  TEBMS  FOE  THE  UTEGEATIOI 

C 


ZL<  D- 

0.0000 

ZL(  2)- 

0.1000 

ZL(  3)- 

0.2000 

ZL(  4)- 

0.3000 

ZL(  5)- 

0.4000 

ZL(  6)- 

O.SOOO 

ZL(  7)- 

0.6000 

ZL(  8)- 

0.7000 

ZL(  9)- 

1.0000 

ZL(10}- 

2.0000 

ZL(ll)- 

3.0000 

ZL(12)- 

4.0000 

ZL(i3)- 

7.7143 

ZI.(14)- 

10.0000 

ZL(IS)- 

16.1429 

ZL(16)- 

16.0000 
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ZL(17)«  17.0000 
XL(18)»  18.0000 
1L(19}-  19.0000 
XL(30)«  20.0000 

n.(21)»  20.1000 

ZL(22)-  20.2000 
XL(23)"  20.3000 
XL(24)«  20.4000 
XL(26)>i  20.4167 

00  102  V  -  1.26 

XL(I)  -  tt/20)*XL(I)  -  L/2 


coirnuE 

BR(  !)• 

0.000 

BR(  2)« 

0.486 

me  3)« 

0.868 

Wl(  4)- 

0.778 

n(  6)« 

0.871 

n(  6)« 

0.946 

Ml(  7)- 

1.010 

BR(  8)- 

1.080 

BRC  9)« 

1.180 

BR(IO)- 

1.410 

BR(11)« 

1.670 

BR(12}« 

1.660 

Ba(13)« 

1.670 

BE(14)- 

1.670 

BR(16)« 

1.670 

BK(18)« 

1.630 

m(17)- 

1.370 

tt(l8)« 

0.919 

n(19)« 

0.448 

BR(20)- 

0.196 

BE(21)« 

0.188 

n(22)- 

0.168 

BR(23}« 

0.132 

BR(24)- 

0.063 

8R(26)- 

0.000 

DO  103  « 

-  1,26 
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XL(N)  -  ZL(M)/L 
Bi(M)  -  Bt(M)/L 
103  CQITIlUB 


C 

pisn  ^n/DBLT 
isn  -piszN 
c 

c  snrounoi  bbqzis 
c 

DO  100  I-1.I8IM 

c 

C  GiLCUUTB  ns  OtAO  FOECB.  IRinAII  US  BBAG  OTEE  TIE  fEIICLE 

C 

00  101  K-1.26 
OCF-H-ZL(K)*q 
VBCl (K)«BE(K)*OCP*ABS(tKr) 
flC2(K)-IE(K)*1ICP*ii8(TICP)*ZL(K) 

101  ooimoE 

GiLL  lEAPCOS.flCl.ZL.IEiVB) 

GiLL  mF(26,TlC2,XL»PXTGI) 

BifE^Z*IEAfE 

PITC1«GDZ*PITCI 

G 

c  coopluo  eqdatiois 

G 

OIIDf-(IT-liqDOT)*(-HEAfE)i-(IIAS8*Z04^»|DOT)*PnGl 

DqDT«(IIASS>ZVDOT)*PlTCl4()US8*ZO-4ll)DOT)*(-lEATE) 

ClDT«(IT-KQD0T)*IIASS*Za*Q**2-(IUSS*Za<^ZqD0T)*IUS8*Za*H*q 

G20T-*(llASS*Zin>0T)«IIA8S*Za*V*q-»(||i8S*ZQ4|IIID0T)*IIASS*Za*q**2 

c 

OH-D«D?/Df 

oq-DqoT/0? 

Gl-GiOf/D? 

G2«G2D?/DT 

G 

HDOT  -All«H‘M12*q«A13*(ZQ8*C0S(llETA) 

*  ♦ZQB*SIE(THETA})<fDV4Cl 

qOOT  -A21*V<»A22*q-»A23«(ZO*G08(TBSTA) 

«  ♦ZaB*SII(TIETA))«0q-»G2 

TBDOT-q 

IF  (I.0T.3)  00  TO  160 


n 


c 

c  niTiAL  msT  OEon  htbsiatioi 

c 

H  -  W  *  OELT*VDaT 
Q  -  Q  mT*QOOT 

TBBTA  -  THETA  +  OBLT*THZOOT 

C 

ri(I)>HDOT 

F2(I)«QD0T 

F3(I)«THEDOT 

C 

C  ADAHS>BASHrOETI  UTEGBATIOI 

C 

150  Pl(4)*All*»-»A12*q4A13*(XGB*C0S(THBTA)4Z6B*SII(TBETA)  )+DU-i>Cl 
P2(4)-A21*IH>A22*q4A23*<XCB*C0S(TBBTA)<i^ZQB*SIl(THETA))‘^DQ4C2 
F3(4)-Q 
C 

W  -H  4(DKLT/24.0)*(55.0*P1(4)>69.0*P1(3)437.0*F1(2} 

•  •9.0«P1(1)) 

q  -q  4(DELT/24.0)*<56.0*P2(4)*69.0*P2<3)437.0*P2(2) 

•  •9.0«P2(1» 

1HETA*T1BTA4(DBLT/24.0)*(S5.0*P3(4)-59.0«P3(3)437.0*P3(2) 

•  •9.0*P3(i)) 

C 

P1(1)-P1(2} 

P1(2)«P1(3) 

Pl(3)-Pi(4) 

P2(1)«P2(2) 

P2(2)-P2(3) 

P3(3)-P2(4) 

P3(i)-P3(2) 

PS(2)-F3(3) 

P3(3)-P3(4) 

C 

TIIIE-I4DELT 

J-J41 

IF  (J.IZ.IPBIT)  00  TO  100 
A1.FA-H 

ALFA-(-ATAI(ALPA) )*180/PI 

WHITE  (20,991)  TDIE,THETA*180/PI.q,ALPA 

J"0 
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c 


100  COMTiro 
STOP 

Ml  POINAT  (4I1S.6) 

BID 

SOnODTIlB  niPd.A.B.OOT) 

c 

C  lORBKICAL  mUBATIOl  lOOTZn  OSHG  m  TIAPBZOIDAL  IDLE 

C 

Dnfuszoi  Ad)  .1(1) 

ll-l-l 

OOT-O.O 

DO  1  I-l.ll 

00n«O.6*<A(I)+A(I<»l))*(B(m)>B(I)) 

ODT  -OOT^OVTl 
1  ooimoB 
BBTOU 
BID 
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